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INTRODUCTION 


A Monte  Carlo  method  has  been  developed  to  solve  the  time-dependent  heat 
conduction  or  diffusion  equation.  The  method  has  been  implemented  in  a versa- 
tile computer  program,  and  applied  to  obtain  the  solution  of  heat  conduction 
problems  in  complex  geometry.  Other  possible  applications  include  particle 
diffusion  and  neutron  slowing  down. 

The  Monte  Carlo  method  is  a generalization  of  the  "floating  sphere"  method 
1 2 

developed  by  Haji-Sheikh  and  by  Muller  . It  is  shown  that  the  solution  can  be 
estimated  by  constructing  a random  walk  based  on  the  selection  of  position  and 
time  using  probability  distributions  based  on  known  Green's  functions.  These 
Green's  functions  satisfy  appropriate  boundary  conditions  on  the  surface  of 
arbitrary  volumes.  They  can  be  obtained  for  a wide  class  of  volumes^.  The  form 
of  the  known  solution  simplifies  if  the  volume  is  wholly  contained  within  a 
homogeneous  medium,  but  Haji-Sheikh's  restriction  of  the  volume  to  a sphere 
is  not  necessary.  Our  approach  is  to  specialize  the  class  of  volumes  to  rec- 
tangular parallelepipeds  of  arbitrary  size.  This  leads  to  an  exact  solution  of 
the  heat  transfer  problem  if  all  the  boundaries  of  the  configuration  are  planar, 
and  to  a solution  with  any  arbitrarily  preset  degree  of  accuracy  if  curved 
boundaries  are  involved. 

The  method  provides  the  possibility  of  solving  time-dependent  heat  conduction 
problems  with  internal  heat  sources  and  a variety  of  boundary  conditions.  The 
current  program  treats  the  various  boundary  conditions,  but  not  internal  heat 
sources.  An  exact  treatment  of  the  linear  heat  conduction  problem  is  obtained. 

If  the  conduction  properties  depend  on  the  local  temperature,  we  linearize  the 
problem  by  breaking  up  the  calculation  into  small  time  steps,  assuming  no 
temperature  dependence  of  the  parameters  during  the  time  step. 


S 


Complex  three-dimensional  configurations  can  be  treated.  The  geometrical 

4 

description  is  of  the  Combinatorial  Geometry  type.  Geometrical  bodies  are  de- 
fined in  terms  of  intersections  of  quadratic  surfaces.  Geometrical  regions  are 
defined  in  terms  of  intersections  of  bodies. 


II.  THEORY 


L*  i 
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1 
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1.  The  Heat  Conduction  Equation 

Let  us  consider  the  heat  conduction  equation: 


‘ p£  ^VT<x,t)  + = ~ Q<x,t) 


(1) 


T(x,t)  is  the  temperature  at  point  x,  time  t;  K is  the  thermal  conductivity, 
P is  the  density,  c the  specific  heat,  and  Q the  heat  source  density.  The 
problem  is  defined  for  xefl  where  Q is  a volume  surrounded  by  a surface  Z. 
K,p,c,  are  continuous  functions  of  position,  except  across  specified  internal 
boundaries.  The  boundary  conditions  on  T are  specified  as  follows: 

- Internal  boundaries: 

T continuous  j 
Kn»VT  continuous  I 


(2) 


External  boundary  Z 

The  external  boundaries  conditions  we  will  consider  are  either: 

(3a) 


T (x, t)  known  for  xel,  0<t<tQ 


or:  i?-KVT(x,t)  = h(TG<x,t)  - T(x,t)) 


h,  T_(x,t)  known 
G 


for  xeZ,  0<t<tQ 


(3b) 


where  rf  is  the  outward  normal  to  the  external  boundary  at  point  x,  and  h is 
the  coefficient  of  surface  heat  transfer. 

- Initial  conditions: 

T (x, 0)  known  for  xefl  (4) 


2.  The  Green's  Function  Method 


. 1 
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The  problem  we  consider  consists  of  finding  the  temperature  T(x0,tQ)  at 


a point  xQ  eft,  to>0* 


Let  us  surround  the  point  x^  by  a volume  Vefi , the  volume  V being  bounded 


by  a surface  S,  and  define  a Green's  function  G(xQ-x,tQ-t)  for  xeV,  t<t0„  It 


satisfies  the  differential  equation 


, 3G(x  -x,t  -t) 

V-kV  — G(x  -x,t  -t)  + — 

pc  o o at 


= o , 


(5) 


the  boundary  conditions  are  as  follows: 

- Internal  boundaries  - or  portions  (if  any)  of  the  internal  boundaries  of  the 


full  configuration  which  are  within  V: 


— G continuous 
pc 


K 7 — G continuous 
pc 


(6) 


- External  boundary  S 


G (Xg-X/tg-t)  = o for  xeS  , 0<t<tr 


- - 0 


(7) 


- Initial  condition 


G (Xq-x,0)  = 5 (x-xQ) 


(8) 


A particular  linear  combination  of  Equations  (1)  and  (5)  can  be  written 
in  the  form 


Jo  °dt  /v  ||T(x't)  VK7  £ G(xo"x'Vt 


+ G(x0-x,tQ-t)  — VKVT(x,t)f  dVx 


(9) 


Jo°dtl 

/ot0dt/vlTU' 


G(X0_X'Vt)  ^’2(X't)  dVx 


t) 


3G(Vx'Vt) 


at 


f G(x0-xft0-t)  f dvx 


'll 


L_ 


The  first  volume  integral  of  the  left  hand  side  can  be  reduced  to  a 
surface  integral  by  applying  a generalized  Green's  theorem.  The  time  integra- 
tion of  the  right  hand  side  can  be  performed.  One  obtains: 


J. 


°dt  | ^-T(x,t)  KV  G(x  -x,t  -t) 

0 Js  I ec  0 0 


+ G(x0-x,tQ-t)  — KVT(x,t)^  *dSx 


+fo  ^ fvG(X o'X'Vt)  & Q(X' 


t)  dV 


(10) 


T(x,tQ)  G(xq-x,0)  dVx 


P 


T(x,0)  G (x0-xftQ)  dVx 


The  second  term  of  the  surface  integral  vanishes  because  of  the  boundary 
condition  (7) . The  first  volume  integral  of  the  right  hand  side  can  be  evaluated 
taking  the  initial  condition  (8)  into  account.  One  therefore  obtains,  after  re- 


arranging the  terms,  and  substituting  t=tQ-T: 


’^O'V  =/Q  °dt  y^G(xQ-x,T)  ^Q(*,t0 
■/. 

V I °dt  / -T (x , t— T ) KV  — 

Jo  Js  0 pc 


•T)  dV 


+ / T (x,0)  G (x0-x,tQ)  dVx 


(11) 


— G (xQ-x, t) *dSx 


Equation  (11)  can  be  considered  as  an  integral  equation  for  T(xo,tQ), 


v~ 


— 


The  choice  of  the  volume  V is  arbitrary,  provided  Ve(2.  If  V=£2  and  the 
corresponding  G is  known,  the  problem  of  finding  T(Xg,tQ)  is  reduced  to  quadra- 
tures (assuming  that  T is  known  on  the  external  boundary  E=S) • In  practice, 
the  choice  of  V is  limited  to  volumes  for  which  the  Green's  function  G is  known, 
or  readily  computable,  such  as  spheres,  rectangular  parallelepipeds,  etc. 

3.  The  Monte  Carlo  Method 

A conceptually  simple  Monte  Carlo  technique  can  be  constructed  for  the 
solution  of  Equation  (11) . Before  describing  it,  a lemma  has  to  be  proven. 

Let  us  consider  an  arbitrary  volume  V,  with  no  internal  heat  sources 
(Q=0) , and  boundary  conditions  T(x,t)=l  for  xeS,  t>0  and  T(x,0)=l  for  xeV. 

The  solution  is  then  T(x,t)=l  for  any  internal  point  x,  at  any  time  t>0,  in 
particular  at  xQ,tg. 

Substituting  this  into  Equation  (11)  we  obtain: 


1 = 


f G(x0-x,to>dvx  ♦ j\,  f-™-±six0-x. T)-  as. 


(12) 


Letting  tQ->-®  in  Equation  (12)  we  get: 


1 = 


"Kv  ^G(xo'x'T),dss 


as  lim  < I G(Xg-x,tQ)dVx^  ->0 
tg-x»|  JV 


Substituting  that  result  into  Equation  (12)  we  obtain: 


-/v 


G(x0-x,t0)  dVx 


(13) 


A Monte  Carlo  algorithm  for  the  solution  of  Equation  (11)  is  then  as 
follows. 

First  estimate  the  time-volume  integral  involving  internal  sources  (if 
any).  For  estimating  the  remaining  two  terms,  first  select  a time  t>0  from 


P(T)dT 


_1 

pc 


G (Xg“ X , X ) 


•3s 


dt 


(14) 
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If  T>t g (which,  according  to  13,  happens  with  probability  J^G  (x0"X't0)dVx)' 
the  second  term  of  Equation  (11)  has  to  be  sampled:  sample  a point  x in  V with 

density  proportional  to  G<xQ-x,t0),  and  score  the  (known)  T(x,0)  as  a contri- 
bution to  T(xQ,tg)  - the  history  terminates.  If  T<tQ,  the  last  term  of  Equation 

1 

(11)  has  to  be  sampled.  Sample  -KV  — G(xQ-x,T)*n  for  a point  x on  S.  (n  is 
the  outward  normal  to  S at  x;  T(x,tQ-T)  is  the  estimate  for  T(xuft0).  Two  cases 
can  occur:  in  one  case,  the  sample  x is  on  a common  part  of  S and  Z,  then 
T(x,tQ-T)  is  either  known  (if  Equation  (3a)  applies)  or  obtainable  (if  Equation 
(3b)  applies)  as  described  in  Section  II. 5 below,  and  the  random  walk  terminates. 
In  the  other  case,  x is  on  S but  internal  to  Z.  The  point  x is  then  surrounded 
by  another  volume  Vtfl,  and  the  procedure  to  estimate  T(x,tQ— t)  is  identical  to 
that  just  described  to  estimate  Ttx^jt^).  In  all  cases  the  random  walk  terminates 
when  a known  temperature  is  encountered  (either  at  t=0  or  on  the  boundary) . 

4.  Inhomogeneous  External  Boundary  Conditions 

The  case  of  inhomogeneous  boundary  conditions  can  be  treated  in  the  same 
fashion  as  the  case  of  known  temperature  condition,  provided  some  adjustments 
are  made. 

If,  at  any  step  of  the  random  walk,  the  surface  S surrounding  V has  a part 
(S=S^  + S^)  in  common  with  a portion  of  Z where  inhomogeneous  boundary  condi- 
tions apply,  the  boundary  conditions  imposed  on  G (Equation  7)  either  have  to, 
or  can  be  modified  to: 

G(x0-x,tQ-t)  = 0 for  xeS^,  0<t<t  (15a) 

it-KV  ^ G(x0-x,t0-t)  = - h ^ G(x0-x,t0-t) 


for  xeS2,  0£t<tQ 


(15b) 


Equations  (9)  and  (10)  are  still  valid.  The  first  term  of  Equation  (10)  is 
rewritten  in  shorthand  notation  as : 


The  second  term  of  the  - integral  vanishes  because  of  the  boundary  condition 
(15a).  In  the  S -integral,  we  replace  T by  its  expression  obtained  from  (3b), 
and  G by  its  expression  obtained  from  (15b) . We  obtain 


F-/ot0dt/s;Tls-a-3s+/ot0at/s2l; 


(T  - rn-KVT  )K  — G 


G h 


- (~-  KV  — G)  — KVT  } *dS 
h pc  pc  ' 


Combining  terms,  one  is  left  with 


F ■=  / °dt(  / -T  KV  — G'dS  + / -T  KV  ~~  G-dS  ’ 

0 )Js1  J S2  G pc 


The  equivalent  of  Equation  (11)  is  then 


pc 


T (xn,t  ) = / °dT  / G — Q dV  + | T (x, 0)  G dV 

0 Jo  Jv  pc  Jv 


"dT<,  I -T  KV  — G*dS  + -T  KV  ~ G-^s! 

pc  lc  G pc 


(16) 
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Equation  (11)  collapses  if  xQ  is  on  S.  Equation  (16)  collapses  if  xQ  is  on 
, but  is  well-behaved  if  xQ  is  on  S2  (provided  h<») . This  shows  that  one 
may  use  Equation  (11)  (with  simpler  Green's  functions)  until  the  necessity 
occurs  to  estimate  the  temperature  at  an  external  boundary  where  inhomogeneous 
boundary  conditions  are  imposed.  At  that  point,  one  has  to  switch  to  the  more 
complex  boundary  condition  (15b) , which  leads  to  Equation  (16) . 

The  Monte  Carlo  procedure  corresponding  to  Equation  (16)  is  to  estimate 
the  internal  heat  source  term  (if  any)  first.  Then  select  t from  Equation  (14) . 
If  x>tg,  sample  the  second  term  and  terminate  the  history.  If  T<tg,  sample  the 
surface  terms.  If  xeS_,  estimate  T_  and  terminate  history.  If  xeS, , estimate 
T(x,t0-T)  either  as  a known  temperature , or  by  estimating  it  using  Equation  (11). 

5.  Linearization  of  a Non-Linear  Conduction  Equation 

The  heat  conduction  equation  (1)  is  non-linear  if  the  thermal  properties 
(K,p,c)  depend  on  the  temperature  T.  Consider  a (T)  which  depends  explicitly 
on  temperature  only. 

Let  us  introduce  the  new  variable 


a (T)dT 

It  follows  from  (17)  that 


(17) 


§f-a£  , W-.VT 


(18) 


Substituting  (18)  into  (1)  we  obtain 


_L  v^70+i||  Q 

pc  a a 3t  pc 


or: 


— V - V0  + ^ = -*  Q 

pc  a dt  pc  * 


(19) 
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If  the  problem  at  hand  is  such  that  one  can  choose  an  "a"  so  that  K/a  is 


i 


r. 


* 

£ 

“ 


i 

t 

f 

> 


time  independent  and  position  independent  in  any  region  (thus  exhibiting  only 
discontinuous  variation  across  specified  boundaries) , Equation  (19)  can  be 
simplified  to 


- ^0  + 
pc 


30 

at 


(20) 


In  practice  this  situation  occurs  in  problems  involving  a single  homogeneous 
material.  If  the  configuration  is  heterogeneous,  we  propose  to  choose  a value 
of  a such  that  K/a  is  constant  for  the  most  important  material.  Regions  involv- 
ing other  materials  can  be  subdivided  into  small  enough  subregions  so  that  K/a 
does  not  vary  appreciably.  In  that  case.  Equation  (20)  holds  as  a reasonable 
approximation. 

In  any  case,  the  diffusivity  K/pc  is  still  temperature  (and  therefore  ©)- de- 
pendent. We  propose  to  subdivide  the  entire  configuration  into  small  enough  re- 
gions so  that  the  diffusivity  K/Pc  can  be  reasonably  approximated  as  being 
constant  within  each  region. 

The  time-dependence  of  the  diffusivity  can  also  be  handled  approximately 
by  splitting  the  time  span  from  t=0  to  t*tQ  into  time  steps  short  enough 
so  that  the  thermal  parameters  do  not  vary  appreciably  over  each  time  step. 
Knowing  the  temperature  at  the  beginning  of  the  time  step,  the  thermal  parameters 
can  be  calculated  and  considered  as  constant  throughout  the  time  step.  The 
temperature  has  to  be  calculated  everywhere  at  the  end  of  the  time  step  - and 
the  new  temperature  distribution  can  serve  as  initial  conditions  for  the  next 
time  step. 

If  the  time  span  has  been  subdivided  into  small  enough  time  bins,  and  the 
homogeneous  regions  have  been  subdivided  into  small  enough  subregions  as  dis- 
cussed above.  Equation  (20)  applies  with: 

K/a  = const 
Pc/a  = const 

in  each  subregion  and  in  each  time  bin. 
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Equation  (20)  is  identical  in  form  to  Equation  (1)  provided  one  makes  the 
following  substitutions : 

0 by  T 

K/a  by  K (21) 

Pc/a  by  pc 

The  only  remaining  problem  is  that  of  treating  the  inhomogeneous  boundary 
condition  (3b).  It  is  known ^ that  the  transformation  (17)  cannot  be  performed 
exactly  in  the  presence  of  such  inhomogeneous  boundary  conditions.  One  can, 
however,  write  the  boundary  condition  in  the  form 

[v-'«  -•<-«]  M) 

where  s 

0G(x,t)  = 0 (x , t ) + a(TG(x,t)  - T (x, t) ) (23) 

<\( 

The  problem  is  that  0G(x,t)  is  not  a known  function,  as  it  involves  the  unknowns 
0 (x , t)  and  T(x,t).  We  propose,  however,  the  following  approximation,  which  is 
in  the  spirit  of  all  those  made  in  this  section: 

0G  (x,t)  = 0 (x,0)  + a (Tg  (x,t)  - T (x,  0) ) (24) 

Equation  (22)  is  then  identical  to  Equation  (3b)  when  substitutions  (21)  to- 
gether with  the  replacement  of 
h/a  by  h 

(25) 

<v 

0 by  Tq 


are  made. 


I 


r;  • 


III.  GREEN'S  FUNCTIONS  IN  HOMOGENEOUS  RECTANGULAR  PARALLELEPIPEDS 

The  Monte  Carlo  method  described  in  Section  II  is  valid  for  any  volume  V 

surrounding  xQ,  provided  V is  entirely  within  the  configuration  unaer  examination. 

In  practice,  the  choice  of  volumes  V is  limited  to  such  shapes  for  which  the 

Green's  functions  are  known  or  easily  computable.  Muller1  considered  a variety 

2 

of  shapes  for  the  solution  of  the  Dirichlet  problem.  Haji-Sheikh  considered 
only  spheres  for  the  solution  of  heat  conduction  problems.  Carslaw  and  Jaeger3 
give  Green's  functions  for  a variety  of  shapes. 

In  this  section,  we  give  all  the  relevant  properties  of  Green's  functions 
defined  over  rectangular  parallelepipeds.  The  restriction  of  volumes  V to  be 
rectangular  parallelepipeds  permits  an  exact  solution  in  the  case  of  configura- 
tions with  piece-wise  planar  boundaries,  or  solutions  to  an  arbitrary  degree  of 
accuracy  if  curved  boundaries  are  involved.  Throughout  this  section,  we  assume 
that  K and  the  combination  pc  are  constant,  and  introduce  the  coefficient 
of  diffusivity 

D = K/pc.  (1) 

1.  Separation  of  Variables 

We  propose  to  solve  for  G(x^,x2,x3»t)  which  satisfies 


n [~8  G + 9 G 3 G ~|  3G  _ 

L 2 * 2 „2_J  " 3t 

^x_L  3x2  3x3 

for  xeV,  where  V is  defined  by 


(2) 


ai  <x^  , i = 1,3 


and  t>0. 
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with  the  boundary  and  initial  conditions 


x2,x3,0) 

= 5(x1)5(x2)6(x3) 

(3) 

x2,x3,t) 

+ 3G(XwX  ,x-,t)  + 

= - 0.“  = — =— , x.  = a.  , i=l, 3 

pi  3x.  '1  1 

1 

(4) 

+ 

, 0~  are 

given  constants: 

+ 

+ aT  >_  o 

i = 1,3 

(5) 

+ 

+ 3’  > 0 

i = 1,3 

(6) 

The  inequalities  (5)  insure  that  the  point  xi=0  is  inside  V.  If  all  0^ 

are  zero,  the  boundary  condition  (4)  is  tailored  to  that  of  Equation  (II. 7). 

+ 

Positive  values  of  0~  cam  be  chosen  to  reproduce  the  boundary  condition  (II. 15b), 
Let  us  assume  that  the  solution  G can  be  written  in  the  form 
G(xlfx2,x3,t)  = X1(x1,t)X2(x2,t)X3(x3,t)  (7) 

Substituting  (7)  into  (2)  we  obtain 


1 r 3 Xi  3Xi  1 

X1X2X3  * r LD  — - arj  = ° (8> 

1=1  i 3x . 

i 

A solution  of  (2)  with  the  conditions  (3),  (4)  is  therefore  (7)  with 


3 X.  3X. 

_ i i 

D - = 0 

3x.2  3t 

i 


i = 1,3 


with  the  initial  amd  boundary  conditions 

Xi(xi,0)  = SfxJ  i = 1,3  (10) 

+ 3X  (x  ,t)  + 

X.(Xi,t)  = -0-  — , x.  = u.  , i = 1,3  (11) 

The  four-dimensional  problems  (2-4)  have  been  reduced  to  three  independent  two- 
dimensional  problems  (9-11) . 
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As  discussed  in  Section  IX. 3,  one  first  needs  to  sample  a time  t from 


the  distribution  (Equation  1.13-14) 
/t“p(t)dt  = y*G(X1,X2,X3,t)  dXjj 


^2^3 


Using  the  expression  (7)  for  G,  we  obtain 


(12) 


r.  r*  *2  [' 

t P(t)dt  = / dxx  /_  X2(x2,t)  dx2  /_ 

J -'a"  J a~ 


2 / - X3<x3»t)  ^3 

a3 


= F1(t)  F2(t)  F3(t) 


(13) 


where 


Fi(t) 


f °+ 

= / _ 1 Xi(xi,t)  dxi  , i = 1,: 
a . 


(14) 


To  sample  the  cumulative  distribution  (13)  we  can  sample  each  of  the  three 
cumulative  distributions  (14)  and  retain  the  smallest  sample. 

Indeed,  the  probability  density  function  (pdf)  of  the  retained  sample  is: 

dF  dF  dF 

pltiat  - — r2r3  . Fl  _ f3  . Fl  _ (15) 

' df  [w3]  SE0- 

As  further  discussed  in  Section  II. 3,  the  sample  t has  to  be  compared  to  a 
given  cutoff  time  tQ.  If  t>tQ,  a position  xlfx2,x3  has  to  be  sampled  from  a 
pdf  proportional  to  Gtx^x^x^tg) . As  G is  separable,  the  problem  is  reduced 


to  sampling  x.  from  X.<x.,t0),  for  i=l,2,3. 


If  the  sample  t is  less  them  tQ,  a point  x^,x2»x3  has  to  be  sampled  on 
S from  a pdf  proportional  to  -Eh^-  g (x3,x2,x3,t) . The  probability  of  such  a 
point  being  on  either  of  the  two  sides  of  the  parallelepiped  which  are  parallel 
to  the  ^2,x3  Plane  is: 


r 3xi  3xi  *i  /~a2  C a 3 

1 = ° L'  ^ + ^ *l'*U  Ja-  VX2't)dx2  Ja-  X3(x3' 


t)dx. 


p 1 32X  (X  ,t)  -| 

• L°  J. 


F2(t)F3(t) 


"l  1 


Using  Equation  (9)  we  obtain 


[ J _ 1 alxi(xi't)dxi] 


F2(t)F3(t) 


dF  (t) 

-dt~  F2(t)F3(t) 


The  probability  P of  selecting  a point  on  either  one  of  the  pair  of  sides  1 is 
therefore  proportional  to  having  selected  the  smallest  time  t from  the  cumula- 
tive distribution  F^(t)  (see  Equation  15).  This  can  easily  be  generalized  to 
the  pairs  of  sides  2 and  3. 

Having  selected  i,  (say  i-1),  x^  is  set  to  or  with  probability 
3X1  - 3X1  + 

- — — ( a, ,t)  or  - — (a  ,t),  respectively.  The  remaining  two  coordinates 

dX_  X dX_  X 


<i=2,3)  have  to  be  sampled  from  X. (x.,t) 


3.  The  One- Dimensional  Green's  Function  Vanishing  at  Boundaries 
We  now  consider  the  case  of  Green's  functions  with  G=0  boundary  conditions 
at  both  boundaries. 


The  Green's  function  satisfies 


3 X(x,t)  3X(x,t)  _ 
.2  3t 

OX 


with  the  boundary  and  initial  conditions 
X (x,0)  = 6 (x) 


X (+a/2,t)  = 0 


therefore  treating  the  symmetric  case  a =-a  =•  a/2. 

Let  us  introduce  the  reduced  variables 
C = x/a 


t = Dt/a 


G(S,f)  = X (x,t) 

Equation  (16-18)  become 


3 G(C.T)  3G(£.t) 


G(C,0)  = 6(5) 
G (+1/2 ,t)  = 0 


The  solution  G can  be  expressed  either  in  the  eigenfunction  expansion: 
G(5,t)  =2  £ cos((2n  + Dirpexp  jj"  (2n  + 1 ) 2 it2  t] , (23) 


or  in  the  image  expansion: 


G (5» T)  = -i—  £ (-1)"  exp  [-<5  + n)2/4t 

2 /ire  n=-®  L-  — 


The  expansion  (23)  converges  rapidly  for  large  t,  whereas  expansion  (24) 


converges  rapidly  for  small  t. 
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by 


4.  Monte  Carlo  Algorithms  for  Sampling  One- Dimensional 
Green's  Functions  Vanishing  at  Boundaries 


The  cumulative  time  distribution  corresponding  to  Equation  (14)  is  given 


L 


1/2 


F(T)  = / G(S,T)d? 

■1/2 


(25) 


Using  the  eigenfunction  expansion  (23)  we  obtain: 

£-(2n+l)2ir2Tj 


P(T)  E ^rexpl-<2n+i>2" 

n=0 


(26) 


corresponding  to  a probability  density  function 


-§7=4*  1 (-l)n(2n+l)  exp  -(2n+l)2TT2T 

n=0  L J 


(27) 


The  expansion  (27)  is  absolutely  convergent  for  t>0,  and  the  absolute  values 
of  the  terms  are  monotonic ally  decreasing  for  t>t(  1 where 


T(1)  = in  (3)^  it  2 = 0.0139 


At  early  times,  an  approximation  is  suggested:  keep  only  the  terms  n=0 

and  n=+l  in  the  expansion  (24).  Substituting  this  expression  into  (25)  one 
obtains 

F(x)  = \ — / e ~ du  (28) 


-i-±  f 

/iT  J 1/ 


'l/4^r 

Compared  to  Equation  (26) , Equation  (28)  gives  at  least  five  place  accuracy  for 
x<0.075. 


To  sample  a time,  one  can  set  a time  breakpoint  (0.014<t  <0.075)  and 


use  the  early  time  approximation  for  t<t  and  the  eigenfunction  expansion  for 
Optimum  computer  times  are  achieved  by  setting  ^=0.05.  At  that  time, 
only  the  first  two  terms  of  the  eigenfunction  expansion  are  non-negligible. 
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i; 

•u 

K 


f 
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The  Monte  Carlo  algorithm  consists  in  precalculating  F(Tl).  A random 
number  £ is  sampled,  if  £>F(t1)»  the  early  time  case  applies.  If  £<F(Tl) 
the  late  time  case  applies. 

4.1  Early  Time  Case 

One  first  has  to  select  t<t.  from  Equation  (28) . Let  u = l/4/T  . 

-u2  2 

Select  u>u^  from  e du  and  set  t=1/(4u)  . 

As  u^  is  of  the  order  of  unity,  an  efficient  technique  to  select  u is  to 

2 

select  5 from  exp  (-2^5)  and  accept  the  sample  with  probability  exp  (-6  ). 
When  £ is  accepted,  u is  set  equal  to  u^S. 

Once  t has  been  selected,  the  coordinate  x has  to  be  sampled. 

If  diffusion  occurs  to  a boundary,  it  occurs  with  equal  probability  to 
either  boundary  as 

f£  (1/2, t)  = - ||  (-1/2, T)  (29) 

If  an  internal  point  £ has  to  be  sampled,  its  density  is  proportional  to 
G(£,t)  as  obtained  from  Equation  (24),  which  can  be  rewritten  in  the  form 


;,x)  = (4tt2t)-1//2  l I exp 

n=0  (_  L" 


(£+2n)V4! 


- exp  (£-2n-l)2/4^J  - exp  (£+2n+l) 2/4^| 

+ exp  £ (£+2n+2)/4^j  (30) 

The  selection  is  performed  as  follows: 

2 

First  n>0  is  sampled  from  exp  (-■ n /4j) . Then  the  largest  integer  n is  found  such 


n = 2n  + £ , £>0 

If  £>1/2,  the  sample  ri  is  rejected. 
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When  n is  accepted,  the  samples  n,  £ are  tentatively  accepted  as  samples 


i 


- 

1 

I 


i 


I 


i 


of  the  first  of  the  four  terms  constituting  the  n-th  term  of  Equation  (30) . 
The  remaining  three  terms  are  taken  into  account  by  a rejection  technique: 


f 

2 2 

} 

= exp  ' 

r 

(5-2n-l)  - (g-2n) 

/4rf 

= exp  j 

r 

» 

|u+2n+l)2  - (£-2n) 

/4t| 

= exp  | 

r 

\ 

f 

(C+2n+2)2  - (5-2n)2] 

/4t| 

The  sample  ? is  accepted  with  probability  l-r^-r^r^  In  practice,  this  is 
broken  up  into  three  steps.  Only  r^  is  calculated,  and  the  sample  is  accepted 
with  probability  1-r^.  If  not  accepted,  and  r3  are  calculated,  and  the  re- 
maining tests  are  performed. 

4.2  Late  Time  Case 

One  first  has  to  select  t<t^  from  Equation  (27).  Keeping  only  the 
first  two  terms  of  the  expansion,  one  obtains: 


dF 

d?  dt  = 4TT 


exp  (-if' 


f)  - 3 exp( 


-9ir2t) 


dx 


Let  x = exp (-x  (t-t^)) 
e = expt-Sx2-^) 
p(x)dx  “ (l-3£x8)dx 


(31) 


or 


p(x)dx  = dx  + | (l-x®)dx 


To  sample  x,  with  probability  (l-3e) / (l-e/3)  one  samples  the  first  term  by  setting 


x equal  to  a random  number.  With  remaining  probability,  (8/3) e/ (l-e/3) , the  second 
term  needs  to  be  sampled,  x can  be  set  equal  to  any  of  nine  random  numbers  except 
for  the  largest  of  these  nine.  This  technique  turns  out  to  be  time  consuming.  A 
more  efficient  technique  consists  in  setting  x to  a random  number,  and  accept  x 

g 

with  probability  1-x  . If  x is  not  accepted,  it  is  mulitplied  by  another  random 
number.  The  product  is  a valid  sample  x. 


23 


Finally,  the  time  variable  is  set  to 


t = t1  - £n  (x)  / tt  . 

Once  t has  been  selected,  the  coordinate  x has  to  be  sampled. 

If  diffusion  occurs  to  a boundary,  it  occurs  with  equal  probability  to 
either  boundary,  as  Equation  (29)  still  applies. 

If  an  internal  point  £ has  to  be  sampled,  its  density  is  proportional  to 
G(£,x)  as  obtained  from  Equation  (23),  keeping  only  the  first  two  terms 


G (£,  T)d£  = 2 jcos  (tt£) exp  (-ir^x)  + cos  (3ir£) exp  (-9ir2x)jd£ 
« cos(ir£)  + e cos(3tt£) 

* cos(tt£)  + e [4cos2  (ir£)  - 3cos  (tt£)1 


where 


e = exp  (-87T  t) 


Substituting  sin(~£)  = 2x-l  into  equation  (32): 
p(x)dx  * (l-3e)dx  + 16ex(l-x)dx 

p (x)dx  = dx  + 6x  (l-x)dx  (33) 

To  sample  equation  (33),  with  probability  (l-3e)/ (l-e/3)  set  x to  a random  number. 
With  remaining  probability,  (8/3) e/ (l-e/3) , sample  three  random  numbers  and  set 
x equal  to  the  one  in  the  middle  in  the  order  of  magnitude. 

5.  The  One-Dimensional  Green's  Function  with  Homogeneous  Boundary  Conditions 
We  seek  the  solution  of  the  equation 

„ _ WjSjtL  . „ (34, 

3x 

with  the  boundary  and  initial  conditions 


X(x,0)  = 5 (x) 


X(0,t)  = S 


9X  (0,  t) 


X(l,t)  = 0 (37) 

therefore  treating  the  case  of  the  source  point  x=0  on  the  boundary  x=0: 
a =0,  a+=l,  with  homogeneous  boundary  condition  at  x*0  (8  >0)  and  zero  boundary 


condition  at  x=a  ($  =0). 


Let  us  again  introduce  the  reduced  variables  (19) : 

5 = x/a 
t = Dt/a2 
G(£,t)  = X (x,t) 


Equations  (34-37)  become: 
32G(5,t)  3g(5,t) 


G(C,0)  = 6 (£) 


3g (0,T ) 


G(0,T)  - Y ■'  (41a) 

G(1,T)  = 0 (41b) 

The  eigenfunction  expansion: 

00 

G(C,t)  = I c sin (a  (5-1))  exp  (-<*  2x)  (42) 

, n n n 

n=l 

satisfies  the  differential  equation  (39)  and  boundary  condition  (41a) , (41b) . 

Substituting  (42)  into  (41a)  we  obtain  the  eigenvalue  equation: 

tan  a = - ya  , n>l  (43) 

n n — 

Taking  (43)  into  account,  one  can  show  that 

A 

JQ  sin(otn(£-l))  sin  (am(£-l))  d£ 

2 

= (1  + y cos  a )/2  for  nv=n 

n (44) 

= 0 for  m^n 

Equation  (44)  shows  that  the  expansion  (42)  is  in  terms  of  an  orthogonal 
set.  A necessary  condition  for  (40)  to  be  satisfied  is 


f0  sin(anU-l))S(Od£  = JQ  sin(an(5-l))  G U,0)d5 


Substituting  (42)  into  (45)  and  taking  (44)  into  account,  one  obtains 
-2  sina 

cn  - ; — r (46) 

1+Ycos  a 

n 

Expression  (42)  is  therefore  the  solution  of  (39)- (41)  provided  (44)  and  (46) 


are  satisfied 


Equation  (48)  can  be  rewritten  as 


. 

! 

i 


( 

I 

» 


: 

i 


r 


l 

> 


or 


tan  nn  = ynir  - yr^ 

rin  = Arc  tan  (ynir  - ynn) 

n If 2f  « . . 


(49) 


(50) 


Equation  (50)  can  be  used  for  an  iterative  method  of  solution,  where  the  m 
iteration  can  be  obtained  from  the  (m-1) ^ using 


th 


m „ . , m-1, 

nn  = Arc tan  (ymr  - ynR  ) 


(51) 


The  iteration  can  be  started  by  using  an  approximate  expression  for 

tan  n ° £ n °/(ir/2  - n °)  and  substituting  in  Equation  (49): 
n n n 


0..  0.  _ „ 0 

nn  - nn  ) = ynw  - ynn 


(52) 


nn  = I Ynir  + yir/2+l  - / (ynir  + yir/2+1)  -2ny  ir 


2 2 /2y 


(53) 


Recapitulating  the  results  of  the  present  section,  the  eigenvalues  an 
can  be  obtained  by  using  the  starting  value  (53)  , iterating  to  convergence 
using  Equation  (51),  and  substituting  the  solution  nn  into  Equation  (47). 

For  y>0  and  n large  enough  an^-(n-l/2)  n- 

5.2  Early  Time  Approximation 

The  expansion  (43)  is  rapidly  converging  for  large  values  of  t-  It 
does  not  converge  for  t=0.  In  order  to  derive  an  approximation  to  G(£,t)  for 
small  values  of  t,  let  us  first  solve  a problem  satisfying  Equations  (39)-(41a), 
and  replacing  Equation  (41b)  by 

G(£,t)  ->-0  as  £-*»  (54) 


') 

;'Hh| 

t. 
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and  denoting  its  solution  G(£,t)  as  Gq(5,t). 

Let  us  define  the  Laplace  transformation  of  GQ. 


F(5,to)  = / e~  Gq  (£  *t  ) dr 


It  satisfies  the  differential  equation 


- UF  = 0 


The  general  solution  of  (56)  satisfying  the  boundary  conditions  (40) 


and  (54)  is 


F(u)  = 


/> 


**  (?+Tl) 


the  function  f(n)  is  to  be  determined  from  the  boundary  condition  (41a) 


-C/uj  f oo  -(£+ti)/io 

+ / dn  f(n>  s—= 

i4o  J 0 Ju 


E'"1*  /. 


dn  f(n)e 


- (£+n)  /J 


for  £ = 0,  or: 


°dn  f (n)e“r'/u=  = 1 2_ 

Y+l/ /a)  1+y/u 


The  left  hand  side  of  Equation  (58)  defines  a Laplace  transform  of  f(n)» 

with  the  transform  variable  s=/^.  The  inverse  transform  of  the  right  hand 

side  is  <5(n)“  ~ e • 

Y 

Therefore 


f(n)  = 5 (n)  - ^ e_Tl/Y 
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Substituting  this  into  Equation  (57)  we  obtain 


Ft&u*  - e-n/y£ 

Y Jo 


The  inverse  transform  of  F with  respect  to  10  is 


dn  e“n/Y  e-(C+n)  /4T 


The  function 

G(C,t)  = Gq(C,t)  - Gq(2-£,t)  (61) 

satisfies  the  differential  equation  (39),  as  well  as  the  boundary  conditions 
(40)  and  (42).  For  small  enough  t,  (2, T) <<Gq (0, t) , and  therefore  from  (61): 
G (0, T)  ^Gq (0, *0  ; for  small  enough  t,  expression  (61)  approximately  satisfies 
boundary  condition  (41) . 


The  function  G(£,t)  defined  by  Equation  (61)  can  therefore  be  considered 
as  an  approximation  to  the  solution  of  Equation  (39) , (41b) . 


6.  Monte  Carlo  Algorithms  for  Sampling  One-Dimensional 

Green's  Functions  with  Homogeneous  Boundary  Conditions 

Starting  from  the  exact  eigenfunction  expansion  (43) , we  obtain  the 

following  expressions  for  distributions  of  conduction  time. 

The  p.d.f.  of  time  passage  through  the  boundary  at  which  G = -y  — is: 

3n 


P0(T) 


. 3G (0, t)  _ 2Vinan  -an  T 

+ = S 2 6 

n 1+  ycos  a 

n 


For  the  boundary  at  which  G=0,  it  is: 


p1(t) 


, 2sina  (-cosa  ) -a  t 
9G(1,t)  _ _ n n n n 

3x  1 ^ 2 e 

n 1 + ycos 


Fq(t) 


2sina 

PQ  (t)dt  = Z 

n a (1+ycos  a ) 
n n 


-a  2t 
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and 


2sina  (-cosa  ) -a  t 
. » - n n n n 

(T)dT  = I - e 

n a (1+ycos  a ) 
n n 


For  t=0,  the  following  equalities  can  be  derived: 


VT)  (1+y) 


F1(t)  = 1-F0(t) 


On  the  other  hand,  starting  from  the  approximate  early  time  expression  of 


G(£,t),  we  obtain  the  following: 


V”  - In'  2 


rf’d^  r 

Ui  ^ Y Jo 


’ -n/v  fa  e"(?+n)  /4t 
dn  e n/Y  / d£e 


F (x)  = -2— 
O'1'  1+y 


1/Y+t/y 


_ r “ 2 

— f e U du 

^ J 1/2/7+  /t/y 


1 + jl  r 00  e-u2 

* Jr.. 


Fi>'’  - In 


Expressions  (68)  and  (69)  have  been  calculated  numerically  for  a range 
of  Yr  and  compared  to  the  exact  eigenvalue  expressions  (64)  and  (65),  keeping 
up  to  900  terms.  At  least  five  place  agreement  occurs  for  x<0.1.  For  that 
value  of  x,  convergence  of  the  eigenfunction  series  occurs  with  only  8 terms 
(all  900  terms  are  required  for  t=0.01).  We  can  therefore  define  a time  break- 
point t1(=0.1)  below  which  the  early  time  approximation  applies,  and  above 
which  the  eigenfunction  expansion  has  to  be  used. 
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We  now  turn  to  describe  specific  algorithms  for  time  and  position 
selection  for  x<x  ^ and  x>x^.  In  both  cases  three  possibilities  have  to  be 
covered:  conduction  to  the  5=0  boundary  at  time  t,  conduction  to  the  5=1 

boundary  at  time  x , and  diffusion  to  an  internal  point  0<5<  1 at  time  x . 

6.1  Case  of  Early  Times 

6.1.1  Conduction  to  the  5=0  boundary 

A time  x has  to  be  sampled  from  FQ(x)  as  given  by  Equation  (68). 

Introducing  the  new  variables 
u = 5/2/t 

v = (5+n)/2/t  (71) 

a)  = 1/2/t, 


Equation  (68)  becomes 


4 rf » -u2 . i f: 

e du di 

/n  y J 0 


ln  e"1^  f " e"v2  dv 

J(  l+n)w 


Ql/y  Too  . 

Multiplying  the  u-integral  by  — — / e d£,  which  is  equal  to  unity,  and 

J 1 

changing  the  variable  n to  ?-l.  Equation  (72)  becomes 

rp.-w-  e'v2av]  m> 

"<  Ji  -J*  Jn  J 


The  probability  density  function  of  is  proportional  to  f (w) 

to  - -i si - f " a5  .-“2  [1  - 5.-'£2-i>“2] 

A y Jl  U 


dF(u) 
dw  : 


Expression  (74)  is  positive  for  all  i;  provided  uj>1.  As  w>a)1  and 
= 1/2/^”  = 1/2/TT  = 1.58,  the  condition  is  satisfied. 
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The  problem  of  selecting  t<t^  from  Equation  (72)  is  therefore  reduced  to 


selecting  from  Equation  (74) , where 


= 1/2>'7t~ 


This  can  be  done  as  follows: 

« r/ y 

1)  select  (;>1  from  e 1 d c; 


2)  select  from  e dui 


3)  accept  the  sample  with  probability  l-^e  ^ w , and 
2 

set  t=1/(2w)  . If  rejection  occurs,  repeat  steps  1,2,3. 


6.1.2  Conduction  to  the  £=  1 Boundary 


A time  t has  to  be  sampled  from  F^t)  as  given  by  Equation  (70). 

The  selection  for  0<t<°°  is  rather  simple: 

2,2..  2 

T = y X / 4v 

-X  2 - 2 

where  p(X)  = e and  p (v)  = — e v 

/F 

However,  what  is  needed  is  selection  in  a finite  range  where  Tl 

can  be  rather  small. 

To  develop  a selection  technique,  let 


= /F/y  , xx  = /r^?y 


Equation  (70)  becomes 


F (x)<*e 


- (u-x)  (u+x) 
e du 


Let  u-x  = v 


F(x)-2  [-  e- 

Jo 


(v+2x) 
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The  probability  density  function  of  x is 


/ * dF 
Px(x)  = " 


✓n  JO 


_ -v (v+2x)  , 
2v  e dv 


(76) 


Let  us  define  the  pdf 


p(*,vU  ^2ve-v(v+2x) 


(77) 


/ir 


p^(x)  is  the  marginal  distribution  of  (77) 


The  marginal  distribution  of  v is: 


P. 


/"x  _ _ 2 r-  -2vx1  “] 

' (v)  = / p (x,v)dx  cc  — e v 1 -e 

Jo  A L J 


(78) 


/it 


If  v is  selected  from  (78) , the  conditional  distribution  of  x is  ratio  of 
(77)  to  (78) : 


• 2vx 

Px  (x}  v)  <c  2v  e 0<x<x1 


(79) 


To  select  v from  (78),  the  following  techniques  are  efficient: 

2 

_ »v 

If  x1  is  not  small  with  respect  to  unity,  select  v from  e , and  accept 

-2v 

with  probability  1-e 

2 


If  x1  is  small,  select  v from  2v  e (v™/x,  P(\)=e  *)*  and  accept  with 
-2v  x. 


* 1 i- 

probability  (1-e  ) /2vx  . The  optimum  breakpoint  value  of  x^  is  /7r/2. 


The  overall  efficiency  of  the  rejection  technique  is  100%  as  x^-^O  and  as  x^-k*; 


it  reaches  a minimum  value  of  54%  at  x^  = Jv/2. 


Once  v is  selected,  the  selection  of  x from  (79)  is  trivial:  select  X 


-X 


from  e , 0<X<2vx^  and  set  x = - X/2v. 


Finally  the  time  is  given  by 


2 2 2,2  2 

T = y x = y \ / 4v  . 


The  probability  is  Ml-vx^).  This  fact  can  be  used  to  avoid,  with  high 
probability,  the  necessity  to  compute  an  exponential. 
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6.1.3  Conduction  to  an  Internal  Point 

The  problem  consists  in  selecting  a position  coordinate  x from  the  pdf 
(61) , at  a given  T. 

This  equation  can  be  written  as: 


G(5 


'T)  '/o' 


-?2/4t  _ -<?+n)2/4T 


‘1 


_ q-(2-C)2/4t+  e-(2-£+n)2/4x| 


(80) 


To  select  5 from  (80)  one  can  first  select  n from  — e , r from  e ^ 


and  accept  the  sample  with  probability 


[l  - e-|iC+n)2-S2]  /4t  - e“  [(2'^V]  /4t  + e-  [(2-?+n)  V]  /4T] 


The  efficiency  of  the  rejection  technique  deteriorates  for  small  values 
of  g = -y/2/7.  In  that  case,  another  algorithm  becomes  efficient. 

2 

With  probability  1/ (1+g/x) , sample  x from  xe  and  u from  2ue  . With 

2 -x2  — -x2 

remaining  probability  sample  x from  x /2  e and  u from  2//^  e . In  both 
cases  set  n = 2/7  gx  and  £ = 2/T  u,  and  accept  the  sample  with  probability 

p'  = p«4x/(  (2£+n)). 


For  optimum  operation  the  first  sampling  technique  should  be  used  for 
g>gQ,  the  second  one  for  g<gQ.  The  optimum  value  of  gQ  can  be  calculated  for 
T«l: 


The  efficiency  of  the  mixed  algorithm  is  100%  as  g+-0  and  It  assumes  a 

minimum  value  of  47%  at  g = gQ. 
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6.2  Case  of  Late  Times 

6.2.1  Conduction  to  Either  Boundary 

I 

The  time  is  sampled  by  sampling  a random  number  r and  solving  the  equation 

rFi(0)  = Fi(T)  i = 0 or  1 (81) 

i = 0 or  1 corresponds  to  boundary  C = or  C = 1,  respectively. 

For  each  sample  (each  value  of  y),  the  eigenvalues  of  coefficients  are  cal- 
culated. The  solution  of  (81)  is  obtained  by  performing  a binary  search:  if 

T1<T<x2r  where  F^(Xj)  is  known  for  j = 1,2,  the  next  approximation  of  t is 
obtained  by  linear  interpolation  between  x^  and  x^.  The  value  (x^)  is  calcu- 
lated numerically,  and  x3  replaces  either  x^  or  x^,  thus  reducing  the  interval 
x1~x2.  The  search  is  terminated  if  F^ (x^)-f^ (x^)  becomes  small  enough  (5%  of  F^fO)). 
If  x2  = 00  (e.g.,  at  the  beginning  of  the  binary  search) , the  "linear  interpolations" 
for  x^  is  replaced  by  selection  from  the  first  term  of  expansion  (3)  or  (4): 

t3  - Ti  - h 109  <r*Fi(v> 

al 

6.2.2  Conduction  to  an  Internal  Point 

At  late  time  x,  the  pdf  of  5 0<£<1  is  proportional  to: 

2sina  -a  2x 

G(C,x)  = Z t— — sin  (a  ?)e  n 

1+ycosa  n 

n n 


The  cumulative  pdf  of  £ is  therefore 


p(e) 


■ /*5  g(c» 

Jo 


2sina  -a  2x 

x)d£  = E -tt— 2-r-  (1-cosa  £)  e n 

s n a (l+7COsa  ) n^ 


(82) 


n 


n 


As  in  the  case  of  time  distribution,  the  selection  is  performed  by  solving 
r-P(l)  = P(£)  (83) 

using  a binary  search  technique. 


L 


| 

r , 

! 

■M 


ri 


j 


Biased  sampling  of  the  Green's  functions  becomes  necessary  if  boundary  con- 
ditions of  the  type 


T - T = - ft  — 

9 3n 

occurs  with  large  values  of  both  T and  r- 

g 

The  natural  probability  of  scoring  is  small,  but  is  much  larger  than 
any  local  temperature  T.  This  leads  to  a high  variance  of  the  estimates.  The  sit- 

'Xi  <\j 

uation  can  be  corrected  by  proper  biasing  if  reasonable  estimates  T^  of  T^  and  T 
of  the  local  tempera ture  T can  be  made. 

Let  us  assume  that  the  RPP  is  oriented  along  the  X^X2X3  axis,  that  zero 
boundary  conditions  are  satisfied  in  the  x^-x^  an<^  X2-X3  P^anes  as  weH  as  on  one 
of  the  x -x_  planes.  The  condition  G = - is  satisfied  on  the  other  X..-X.- 

X A o H i.  z 

plane.  The  procedure  outlined  in  Section  3 is  altered:  the  time  t^  of  diffusion 

to  the  xi-x3  planes  and  to  the  x^-x^  Planes  are  sampled  first.  The  smallest  of 

t , t. , t„  is  determined  and  set  as  t , where  tn  is  the  time  cutoff,  t is  scaled 
0 12  m 0 m 

2 

to  the  frame  of  x,:  x =Dt  /a,  . If  F„(t)  is  the  function  defined  by  Equation  (64) 

3 m m 3 0 

or  (69)  of  Section  6,  then,  with  probability 

* Fo(0|-VV 

the  score  will  be  T . With  remaining  probability  (1-p)  it  will  be  a local  tempera- 

g 

ture:  with  probability 

pa  - V0|-FiV 

the  score  will  be  a temperature  T at  a point  on  the  opposite  side  of  the  x 

cl  w 

interval,  and,  with  probability 


Pi  = 1-P0(T»)  " W (pg  + Pa  + pi  = L) 


the  score  will  be  a temperature  T^  at  a point  internal  to  the  interval  of  x^. 
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Assuming  that  T is  reasonably  close  to  T , and  that  both  T and  T.  are 
9 9 a 1 

reasonably  close  to  T,  reasonable  biased  probabilities  and  associated  weight 
factors  can  be  defined  as: 


Pg  = PgVN 


W 


0, 

N/T 


f\j  r\j 

p = p T/N 

CL  aL 


W = N/T 


'll  T/ 

Pi  = PiT/N 


W.  = N/T 

l 


'Xi 

where  N = p T + (1-p  )T. 

9 9 9 

The  associated  scores  will  be  W T , W T , W.T.,  which  are  all  approximately 

ggaaii  •* 

equal  to  N. 
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IV.  GREEN'S  FUNCTIONS  IN  INHOMOGENEOUS  RECTANGULAR  PARALLELEPIPEDS 

The  Green's  functions  defined  for  homogeneous  RPP  in  Section  III  can  be 
used  provided  the  RPP  is  fully  contained  in  a homogeneous  region.  Another  type 
of  Green's  functions  has  to  be  introduced  when  a point  xQ  is  placed  on  the 
interface  between  two  dissimilar  materials.  As  in  the  previous  section,  we  re- 
strict the  shape  of  the  volume  V to  a rectangular  parallelepiped.  The  inhomogeneity 
is  restricted  to  two  media  with  a planar  interface,  the  interface  being  parallel 
to  a side  of  the  RPP. 

To  be  more  specific,  the  volume  V is  defined  by: 

a . < x.  < a?  i = 1, 3 (1) 

i — i—i 


The  diffusion  parameters  are: 
"1 


K = K, 


P„  = 
c c. 


K = K„ 


P„  = P„ 
c c„ 


for  x3>0 


for  x3<0 


(2) 


Equations  (5)- (8)  of  Section  1.2  can  be  rewritten  as: 

_K  / 32G  32G  32G  | 3G 

pc  . 2 . 2 , 2 " 3t 

\ 3X1  3X2  3x3  / 

with  the  boundary  and  initial  conditions 
G(x1,x2,x3,0)  = 6 (x1) <5 (x2) S (x3) 
G(x1,x2,x3,t)  = 0 , xi  = 


1 

PC! 


1 3 

pc3  3x3 


pc. 


1*2' 


(3) 


pc. 


X35 

(4) 

, i - 1,3 

(5) 

,x2,0  ,t) 

(6) 

G(xlfx2,0  ,t) 

(7) 
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Unfortunately,  separation  of  variables  does  not  apply  to  this  problem-  Par- 
tial separation  can  be  achieved  in  cylindrical  symmetry  (if  the  RPP  is  replaced 
by  a cylinder  around  the  x^-axis) . The  remaining  problem  is  three-dimensional 
(r,z,t-dependent) , and  an  eigenfunction  expansion  can  be  written  down.  Attempts 
to  implement  an  approach  based  on  this  eigenfunction  expansion  have  been  given  up. 
The  solution  of  the  eigenvalue  problem  is  too  time  consuming  to  be  carried  out 

during  the  course  of  a Monte  Carlo  calculation.  Precalculation  of  sampling  tables 

+ 

for  a representative  set  of  parameters  (K^/K^pc^pc^  cu)  appear  to  be  impractical 
as  the  tables  would  be  too  bulky  to  keep  in  computer  core. 

An  approximation  (which  can  be  reduced  to  any  degree)  has  been  introduced. 

It  is  well  known  that  the  diffusion  equation  can  be  considered  as  an  approximation 
of  the  transport  equation  in  the  limit  of  small  mean  free  paths.  Conversely,  the 
diffusion  equation  can  be  approximated  by  a transport  equation  with  small  mean 
free  paths. 

Consider  the  transport  of  radiation  in  a homogeneous  medium  with  the  following 

properties:  The  root  mean  square  free  path  of  a particle  is  e,  the  velocity  is  vr 

scattering  is  isotropic.  For  small  enough  e,  the  diffusion  approximation  applies, 

* 

with  a diffusivity 

D = K/pc  = ve/6.  (8) 

The  distribution  of  free  paths  is  irrelevant  in  the  limit  of  small  e.  We  are  free 
to  take  constant  free  paths  of  length  e. 

At  the  interface  between  two  different  materials,  the  scattering  is  anisotropic 
with  probability  p , it  is  isotropic  into  medium  1,  and  with  probability  1-p^,  it 
is  isotropic  into  medium  2.  In  order  to  satisfy  the  boundary  conditions  (6) , (7) 


/KjTpc^ 


/ ("V 


1*PC1 


+ /K 


2’P  c2> 


(9) 


In  neutron  transport,  the  root  mean  square  free  path  e = 2X,  where  X 
is  the  mean  free  path,  giving  the  familiar  D = vX/3. 


The  accuracy  of  the  method  is  determined  by  the  smallness  of  e.  Once 


e is  chosen  for  medium  i,  i=l,2,  the  velocity  v.  is  determined  from  Equation  (8). 
i 1 

t l 

An  adequate  accuracy  is  achieved  by  setting  e.  to  the  smallest  of  e.  , z^  , 
where 

•i*  - /6  D. t0/25 

E = L./6  where  L.  is  the  smallest  of 
i a/  i 

t t + 

+a^  i ±<*2  ' a3  ^ ~a2  ^ l-2. 


Setting  e . - will  insure  25  collisions  before  reaching  a given  time  cutoff 
t . Setting  z-  = z-^  will  insure  a number  of  collisions  of  the  order  of  25  before 


leaking  out  of  the  RPP. 


Once  the  particle  diffuses  away  from  the  interface,  one  can  switch  to  the 
floating  homogeneous  RPP  method.  Optimum  computing  times  are  achieved  if  this 

z 

switch  occurs  when  the  particle  diffuses  a distance  of  from  the  interface. 


V. 


HEAT  TRANSFER  GEOMETRY  PACKAGE 


1.  Geometrical  Description 

4 

The  geometrical  description  is  of  the  combinatorial  type  . Regions  are 
described  in  terms  of  intersection  of  bodies.  Bodies  are  described  in  terms 
of  intersection  of  quadratic  surfaces.  The  only  surfaces  currently  being  im- 
plemented consist  of  planes,  spheres,  cylinders,  and  circular  cones  (corres- 
ponding to  combinatorial  bodies  RPP,  BOX,  SPH,  RCC,  TRC,  WED,  ARB,  but  exclud- 
ing REC  and  ELL) . 

If  a configuration  is  axially  repetitive,  only  a single  repetitive  element 
needs  to  be  described,  as  explained  in  Sect .o:  -TI.3.  The  geometrical  package 

is  completely  general.  However,  if  the  cni  . ; iration  exhibits  at  least  partial 
cyclindrical  symmetry,  more  efficient  Teat  transfer  calculations  will  be  per- 
formed if  the  eutis  is  in  the  z-direction  (only  subroutine  SOUSET  is  affected. 
See  Section  VII. 5). 

2.  Geometrical  Input 

The  input  consists  of 

1)  Title  Card 

2)  Surface  Data 

3)  Body  definition  in  terms  of  above  surfaces 

4)  Region  definition  in  terms  of  above  bodies 

It  is  envisioned  to  develop  a preprocessor  which  will  accept  as  input  com- 
binatorial bodies. 

2.1  Title  Card  (I4,19A4) 

IPR  If  zero  - input  and  processed  data  will  be  printed  back. 

HOLL  Any  hollerith  information  serving  as  a title. 

2.2  Surface  Input  (2X,  A3,  8E9.3) 

The  input  number  e required  for  each  surface  can  be  regarded  as  a surface 
thickness.  Its  meaning  is  described  in  Section  V.4.1. 


Plane: 


PLN  c V]_  V2  V3  Hl  H2 


H 
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where  V is  a point  on  the  plane,  and  H is  orthogonal 
to  the  plane,  pointing  "inside"  the  plane  (arbitrary 
normalization) 

SPH  e Vx  V2  V3  R 

where  V is  the  center  and  R the  radius 
CYL  e V;L  V2  V3  Hx  H2  H3  R 

where  V is  a point  on  the  axis,  H points  along  the 
axis  (arbitrary  normalization)  and  R is  the  radius 
CON  e V1  V2  V3  Hx  H2  H3  T 

where  V is  the  vertex,  H points  along  the  axis 
(arbitrary  normalization) , and  T is  the  tangent  of 
the  angle  of  aperture 
End  of  surface  input  marker:  END 

2.3  Body  Input  (2X,  A3,  1515) 

GEN  N N2  N3  ... 

where  |n. | is  surface  number  N^>0  implies  that  the  body  is  inside 

the  surface,  N,<0  - outside.  The  first  N.  = 0 implies  the  end  of 
i r 

list  at  i-1.  If  none,  end  occurs  at  i = 16. 

End  of  body  input  marker:  END 

2.4  Region  Input  (2X,  A3,  1515) 

ANY  Nl  N2  N3  ... 

where  ANY  is  any  three  characters  (^  END)  and  jN^]  = body  number. 
N^ >0  implies  the  region  is  inside  body,  N^<0  - outside.  N^  = 0 
end  of  list  marker. 

End  of  region  input  marker:  END 


Sphere : 


Cylinder: 


Cone: 


3.  The  Input  Processing  Routines 


GENI,  GENIP,  GENICK 

These  subroutines  read  card  input  and  store  processed  data  into  a single 
dimensional  array  called  AST  in  floating  point  mode  and  MAS  in  integer  mode. 

The  layout  of  the  information  is  shown  in  Table  1. 

The  surface  data  is  processed  first,  and  stored  sequentially  in  the  first 
available  location  IS  of  the  array.  One  word  is  left  blank  for  future  use.  The 
TYPE  is  changed  to  a numerical  code  as  shown  in  Table  1.  The  vector  H is  normal- 
ized to  unity.  An  exhaustive  search  of  surface  duplicates  is  performed.  Each 
time  a surface  is  read  in,  the  previous  list  is  examined  for  am  identical  surface 
If  a duplicate  is  found,  the  new  surface  is  ignored,  and  the  parameter  e of  the 
old  surface  is  changed  to  the  smallest  of  the  old  and  new  e.  The  address  JS 
constitute  the  only  identifier  of  the  surface.  A temporary  dictionary 
ISLOC(IS)=JS  is  generated,  which  defines  the  surface  identifier  JS  corresponding 
to  the  ordinal  surface  number  IS.  In  the  case  of  a new  plane  identical  to  an 
old  plane  except  for  the  vector  H pointing  in  the  opposite  direction,  the  old  H 
is  left  unchamged,  but  the  dictionary  entry  of  the  new  plane  is  tagged  with  a 
negative  sign:  ISLOC (IS) =-JS. 

The  body  definition  is  read  in  next.  The  information  is  stored  in  terms 
of  surface  identifiers  as  determined  from  the  surface  dictionary  ISLOC.  A tem- 
porary body  dictionary  IBLOC  is  generated,  which  defines  the  body  identifier  in 
terms  of  body  input  number. 

The  last  block  of  information  read  in  is  the  region  definition.  The  in- 
formation is  stored  in  terms  of  body  identifiers  as  determined  from  the  diction- 
ary IBLOC.  A permanent  region  dictionary  ISLOC  is  generated,  which  defines  the 
region  body  list  identifers  in  terms  of  the  ordinal  region  numbers.  (The  ori- 
ginal surface  dictionary  ISLOC  is  destroyed.) 
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LOCATION 


IB+Jl-1 


TABLE  1 


Layout  of  the  MAS  - AST  Array 


CONTENT 


BLANK 
ITYPE  \ 

V1  J 


BLANK* 

ITYPE 


BLANK 


surface  "1" 


1 = SPH  2 = CYL  3 = CON  4 = PLN 


-cosct 


Body  "IB” 
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TABLE  1 (continued) 


LOCATION 
IPB  = ISLOC(IR) 
IRB+1 
IRB+2 

IRB+m-1 

IRB+m 


CONTENT 


BLANK 


iIBi 


iIB2 


V 


+ IB 
- m 


region  IR 

body  li3t 


IRS  =■  IBLOC(IR) 
IRS+1 
IRS+2 
IRS+3 
IRS+4 
IRS+5 
IRS+6 
IRS+7 


region  IR 
surface  list 


ISM  = IRS+4 (n-1) 
ISM+1 
ISM+2 
ISM+3 
IRS+4n 


JS 

n 

BLANK 

BLANK 

BLANK 

° J 


* 

All  blank  locations  are  zero'd  in  at  input  time 


r = - 1 

The  final  task  of  the  preprocessing  routines  is  to  build,  for  each  region, 
the  list  of  surfaces  involved  in  defining  the  region  boundary.  The  list  is  in 
terms  of  surface  identifiers,  and  is  stored  in  the  MAS  array.  A permanent  re- 
gion dictionary  IBLOC  is  generated,  which  defines  the  region  surface  list  iden- 
tifiers in  terms  of  the  ordinal  region  numbers.  (The  original  body  dictionary 
* [3 

IBLOC  is  destroyed.) 

Throughout  the  input  processing,  a rather  thorough  search  of  errors  and  in- 
consistencies is  performed.  Error  or  warning  messages  are  issued.  If  possible, 
fatal  errors  are  temporarily  fixed  up  to  allow  continuation  of  scanning.  If 
any  fatal  errors  occurred,  a stop  is  executed  at  the  end  of  input  processing. 

4.  The  Geometry  Subroutines 

4.1  The  Box  Fitting  Subroutine  GBOX 

The  principal  purpose  of  this  subroutine  is,  given  a point  in  a region,  to 
construct  a box  centered  around  the  point,  and  fully  contained  within  the  region. 

« 

The  secondary  purpose  is,  given  a point  on  the  boundary  of  a region,  to  construct 
a box  fully  contained  within  the  region,  with  one  face  centered  around  the 
point.  In  both  cases  the  box  has  to  be  "as  large  as  conveniently  possible". 

The  ill” defined  statement  in  quotation  marks  can  be  qualified  as  follows.  The 
efficiency  of  the  diffusion  code  requires  on  one  hand  that  the  smallest  dimension 
of  the  box  be  as  large  as  possible,  and,  on  the  other  hand,  that  the  largest 
fraction  of  the  surface  area  of  the  box  coincides  with  the  region  boundary.  The 
efficiency  of  the  box  construction  routine,  however,  requires  avoiding  lengthy 
calculations  and  tests.  The  subroutine  GBOX  implements  a reasonable  compromise 
between  these  conflicting  requirements. 

The  requirement  that  a fraction  of  the  surface  area  of  the  box  coincides, 
at  least  occasionally,  with  the  region  boundary  is  a necessary  one  for  termina- 
tion of  a Monte  Carlo  history.  This  can  be  achieved  exactly  if  and  only  if  the 
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region  boundary  consists  of  planar  facets.  To  avoid  this  restriction,  we  con- 
sider the  curved  surfaces  as  slightly  diffuse.  A point  within  a distance  e 
from  the  surface  is  considered  to  be  on  the  surface,  where  e is  the  input 
number  defined  in  Section  2.2,  which  should  depend  on  the  diffusion  properties 
of  the  materials  on  either  sides  of  the  surfaces. 

The  construction  of  a box  abutting  a single  surface  is  shown  in  Figure  1 
and  2 for  convex  and  concave  surfaces.  The  distance  is  the  distance  of 
closest  approach,  plus  or  minus  e depending  on  the  convexity  of  the  surface. 

is  the  length  of  the  diagonal  of  the  largest  box,  one  face  of  which  cam  be 
considered  as  being  entirely  "on"  the  surface.  R is  the  length  of  the  diagonal 

of  the  abutting  rectangle.  If  R is  the  small  radius  of  curvature  of  the  sur- 

s 

2 2 

face,  R = 4Rge  ~ e • or  R ^ 2</Rg£.  (Rg  = R for  a sphere  or  a cylinder,  Rg  = 00 
for  a plane,  and  Rg  is  determined  as  a function  of  position  for  the  cone). 

The  subroutine  examines  each  surface  of  the  list  IRS.  It  calculates  the 
distance  of  closest  approach,  and,  if  this  turns  out  to  be  the  currently  short- 
est  one,  calculates  also  R,  D^,  as  well  as  Q where  Q is  the  vector  front  the 
point  of  closest  approach  to  the  given  point  X.  The  distance  of  closest 
approach  to  the  next  nearest  boundary  is  stored  as  D3> 

When  the  list  of  surfaces  has  been  exhausted,  a box  is  tentatively  con- 
structed. The  diagonal  of  the  box  is  D^,  the  shortest  of  D2  and  D^.  The  di- 

-► 

mension  of  the  box  along  Q is  2D^.  The  other  two  dimensions  are  equal  to  2D^, 

/ 5 2 

with  D2  * / (D|j  - D,  )/2.  The  address  ISM  (see  Table  1)  for  the  closest  surface 
is  saved. 

The  box  just  generated  is  perfectly  valid.  For  efficiency  purpose,  two 
tests  are  performed  to  investigate  the  possibility  of  a "better"  box: 


j 
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1)  The  box  does  abut  a surface  which  is  part  of  the  definition  of  the 
region  boundary.  But  does  it  abut  that  part  of  the  surface  which  really 
defines  the  boundary?  The  situation  is  described  in  Figure  3.  The  point  of 
closest  approach  XQ  is  tested  using  subroutine  REGCK  described  in  Section  4.2.2. 

If  the  point  XQ  turns  out  not  to  be  on  the  region  boundary,  the  box  is  changed 
to  a (larger)  cube  with  diagonal  D3:  Dj*D2'  =D3//3.  The  vector  Q is  changed 
(arbitrarily)  to  0,0,1.  ISM  is  set  to  0.  However,  if  XQ  is  on  the  real  boundary, 
the  following  test  is  performed: 

2)  Is  the  abutting  side  of  the  box  much  further  than  the  non-abutting 
sides  of  the  box,  which  we  implement  as  D3>2D2?  If  yes,  diffusion  to  the  abutting 
side  is  rather  unlikely,  and  a cube  of  diagonal  D3  is  more  efficient. 
(D3=D2=D3//3,Q=0,0,1,ISM=0) . The  situation  is  sketched  in  Figure  4. 

The  above  discussion  describes  the  construction  of  a box  centered  around  an 
internal  point  X.  To  obtain  that  case,  the  variable  JSM  in  C0MM0N/GE0M/  has  to 
be  <0. 

For  the  case  of  X "on"  (in  the  diffuse  sense)  a boundary,  the  variable  JSM 
has  to  be  set  to  the  identifier  of  the  particular  surface.  The  box  construction 
proceeds  in  much  the  same  way  as  previously  described,  except  that  abutting  is 
forced  to  surface  JSM.  The  tests  to  replace  the  box  by  a cube  are  bypassed. 

The  vector  Q is  ill-defined.  Another  vector,  QIN  is  generated,  which  is  the  in- 
side normal  of  surface  JSM  at  X. 

4.2  The  Point  Checking  Routines 

The  package  consists  of  subroutine  REGCK,  SURFCK. 

The  variables  of  interest  are  found  in  COMMON /GEOM/  and  consist  of 

X(3)  point  being  checked 

IR  region  number 

IRP  neighbor  region  number 

JSM  identifier  of  surface  "on"  which  the  point  is  assumed  to  be 

LOOP  the  G uber  number 


The  Guber  number  should  be  familiar  to  all  combinatorial  geometricians 


It  is  initialized  to  zero  by  subroutine  GENI.  It  should  be  incremented  by 
unity  each  time  X changes  between  calls  to  any  point  checking  routines.  Its 
function  will  become  apparent  in  the  following  subsections. 

4.2.1  The  Surface  Checking  Routine  SURFCKfJS, ILOOP) 

JS  is  the  surface  identifier.  ILOOP  should  be  numerically  equal  to  LOOP. 

(Do  not  use  LOOP  itself  as  an  argument!)  The  subroutine  performs  simple  alge- 

j 

braic  tests  and  determines  whether  X is  "inside"  or  "outside"  the  surface. 

t 

The  output  ILOOP=LOOP  if  X is  inside  JS,  ILOOP=-LOOP  if  X is  outside  JS. 
The  output  ILOOP  is  also  stored  in  MAS (JS)  (labeled  as  'blank'  in  Table  1). 

4.2.2  The  Region  Checking  Routine  REGCK(JRB, JLOOP) 

The  purpose  of  this  routine  is  to  check  if  a point  X is  in  a region  IR 
(JRB=ISLOC (IR) ) . JLOOP  should  be  numerically  equal  to  LOOP.  The  subroutine 
handles  differently  the  case  of  a point  inside  a region  or  "on"  (in  the  diffuse 
sense)  a boundary. 

4. 2. 2.1  Testing  a Point  Inside  a Region 

The  "inside"  case  will  be  discussed  first.  If  the  point  is  not  intention- 
ally 'on"  a boundary,  one  should  set  JSM<0. 

The  output  of  REGCK  is  JLOOP=LOOP  if  X is  inside  IRB,  JLOOP=-LOOP  if  X is 
outside  IRB.  The  output  JLOOP  is  also  stored  in  MAS (IRB). 

The  subroutine  functions  as  follows: 

It  accesses  the  region  definition  in  terms  of  bodies,  and  examines  each 

) 

body  IB  in  turn. 

It  first  examines  MAS (IB)  and  compares  its  absolute  value  with  LOOP.  If 
they  match,  the  body  has  been  tested  before  for  the  same  point  X and  the 
following  set  of  tests  can  be  bypassed. 
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If  | MAS (IB) | does  not  match  LOOP,  the  body  definition  in  terms  of  sur- 


•• 

I 

; 
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faces  is  accessed,  and  each  surface  JS  is  examined  in  turn. 

If  |MAS(JS) | = LOOP  the  surface  has  been  examined  before.  If  not,  a call 
to  SURFCK  is  performed.  In  both  cases  the  sign  of  MAS(JS)  indicates  whether  X 
is  inside  or  outside  JS.  The  examination  continues  as  long  as  the  result  of 
the  tests  match  the  body  definition.  The  first  mismatch  indicates  that  X is  out- 
side the  body;  KLOOP  is  set  equal  to  -LOOP.  No  mismatch  until  the  end  of  in- 
formation (JS=0)  marker  is  reached  indicates  that  the  point  is  inside  the  body; 
KLOOP  is  set  equal  to  +LOOP.  In  both  cases  KLOOP  is  stored  into  MAS (IB). 

Whether  determined  previously  or  by  the  tests  described  above,  the  sign 
of  MAS (IB)  indicates  whether  X is  inside  or  outside  IB.  The  examination  continues 
as  long  as  the  results  of  the  tests  match  those  of  the  region  definition.  The 
first  mismatch  indicates  that  X is  outside  the  region;  JLOOP  is  set  equal  to 
-LOOP.  No  mismatch  until  the  end  of  information  (IB=0)  marker  is  reached  in- 
dicates that  the  point  is  inside  the  region;  JLOOP  is  set  equal  to  +LOOP.  In 
both  cases  JLOOP  is  stored  in  MAS(IRB). 

4.2.3  Testing  a Point  "On"  a Surface 

In  addition  to  X and  LOOP,  the  following  variables  have  to  be  set: 

JSM  identifier  of  surface  "on"  which  X lies 

IR  known  region  number  on  one  side  of  the  surface 

IRP  neighbor  candidate  being  tested  on  the  other  side 

of  the  surface 


IRB  (argument)  IRB=MAS ( IRP ) 


The  body  checking  proceeds  as  usual  unless  surface  JS=JSM  is  encountered. 
If  it  is,  the  testing  of  that  particular  surface  is  bypassed,  and  the  remaining 
surfaces  are  considered.  If  the  final  result  of  the  body  test  indicates  that 
the  point  is  outside  the  body,  either  assumption  of  X inside  or  outside  JS=JSM 
will  not  change  that  result  and  the  logic  of  the  region  checking  is  unaffected. 
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If,  however,  the  result  indicates  that  the  point  is  inside  the  body,  the  point 
X is  actually  on  the  boundary  of  the  body.  If  the  region  IRP  being  checked 
is  equal  to  IR,  the  result  can  only  be  false,  and  is  therefore  set  as  such.  If, 
however,  IRP  ^ IR,  the  body  test  is  tentatively  assumed  to  match  the  region  de- 
finition of  IRP,  the  code  proceeds  to  check  the  next  body,  and  the  logic  pro- 
ceeds along  the  lines  of  Section  4. 2. 2.1.  The  tentative  assumption  involves  a 
definite  assumption  on  the  position  of  IRP  with  respect  to  surface  JSM.  That 
assumption  is  recorded  as  a signed  integer  K.  In  the  event  that,  for  another 
body,  the  opposite  assumption  has  to  be  made  the  neighbor  IRP  is  declared  as 
invalid. 

For  a successful  IRP,  the  subroutine  also  generates  a vector  QE  which 
points  from  IR  to  IRP.  This  is  done  by  setting  QB=+QIN,  where  QIN  is  a vector 
(generated  by  GBOX,  see  Section  4.1)  pointing  into  surface  JSM,  and  the  sign 
is  determined  by  the  sign  of  K. 

4.3  The  Neighbor  Finding  Subroutine  NEBFND 

The  purpose  of  this  subroutine  is  to  find  the  neighbor  IRP  of  a region  IR 
at  a point  X "on"  a particular  surface.  The  surface  is  defined  by  the  variable 
ISM  (set  by  subroutine  GBOX).  Referring  to  Table  1,  the  n-th  surface  of  the 
IRS  list  is  defined  by  ISM=IRS+4 (n-1) . 

For  efficiency  purposes,  the  subroutine  builds  up  a list  of  neighbors  pre- 
viously determined,  for  each  region  across  each  surface.  The  list  is  stored 
in  the  MAS  array.  The  first  entry  (if  any)  is  stored  in  MASCISM+l).  The 
second  one  (if  any)  is  stored  in  MAS(ISM+2).  If  a third  entry  becomes  necessary, 
the  first  available  location  LS  in  the  MAS  array  is  determined,  and  the  entry 
is  stored  there;  the  address  of  LS  of  that  third  entry  is  stored  in  MAS(ISM+3). 

LS  is  set  to  LS+3  to  reserve  a location  for  a possible  fourth  entry,  and  for  the 
location  of  a possible  fifth  entry,  etc.  The  buildup  of  tables  stops  if  the 
MAS  array  has  been  filled. 
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The  detailed  operation  of  the  subroutine  is  as  follows : it  examines  the 


first  candidate  (IRP=MAS (ISM+1) ) . If  IRP=0,  no  candidate  exists  and  the 
following  step  is  bypassed.  If  IRP>0,  REGCK  is  called  for  checking  IRP.  If 
IRP  is  found  valid,  the  subroutine  returns  to  the  calling  program.  If  not, 
the  next  candidate  is  examined,  until  either  a validity  check  or  IRP=0  is  en- 
countered. In  the  latter  case,  a loop  over  all  IRPs  is  initiated.  The  first 
check  performed  is  whether  the  region  has  been  tested  before  (IRB=ISLOC (JRP) , 
|MAS(IRB) | = LOOP?) . If  yes,  the  next  IRP  is  considered.  If  IRP  has  not  been 
tested  before,  the  next  test  consists  in  determining  whether  the  boundary  of 
IRP  includes  the  particular  surface.  If  not,  the  next  IRP  is  considered.  Other' 
wise,  IRP  is  checked  further  by  calling  REGCK.  If  the  region  is  found  invalid, 
the  next  IRP  is  considered.  The  first  valid  IRP  terminates  the  loop.  It  is 
stored  as  a next  candidate.  Before  returning  to  the  calling  program,  the 
fact  that  IR  is  a neighbor  of  IRP  across  the  same  surface  is  also  recorded. 

4.4  The  Region  Finding  Subroutine  REGFND 

The  subroutine  determines  the  region  number  IR  for  a point  X.  Error 
messages  are  printed  if  the  region  is  undefined  or  multiply  defined.  Subroutine 
REGCK  is  called  for  all  regions. 

4.5  The  Tracking  Routines  Gl,  G1P 


Given  a point  X,  a unit  vector  Q,  and 
a region  IR,  the  subroutine  calculates 
the  number  NI  of  intersections,  and 
the  intersections  DI(I),  1=1,  NI  of 
region  IR  and  of  the  line  passing 
through  X in  the  direction  Q. 


The  operation  of  the  subroutine  is  as  follows: 

First,  all  the  intersections  with  all  the  surfaces  which  are  involved  in 
the  definition  of  IR  are  calculated.  These  include  all  intersections  with  IR, 
and  may  include  some  intersections  which  are  not  on  the  surface  of  IR.  Sub- 
routine G1P  first  orders  all  the  intersections,  then  weeds  out  the  extraneous 
ones.  The  point  X is  moved  to  below  the  first  intersection  (thus  to  a point 
definitely  outside  IR) . Subroutine  REGCK  is  called,  with  a special  switch 
JSf^-33333  which  forces  REGCK  to  check  all  surfaces  without  returning  at  the 
first  mismatch  encountered.  A loop  is  initiated,  which  runs  through  all  the 
tentative  intersections.  For  each  intersection,  the  logical  information  of 
whether  X is  in  or  out  of  the  corresponding  surface  is  negated.  Subroutine  REGCK, 
which  queries  that  information,  and  its  logical  output  (in  or  out  of  IR)  is  com- 
pared to  its  previous  logical  output.  If  the  two  match,  the  corresponding 
intersection  is  not  on  the  surface  of  IR,  and  is  weeded  out.  If  the  two  do  not 
match,  a true  intersection  has  been  encountered  and  is  therefore  retained. 
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VI.  THE  MAJOR  MONTE  CARLO  HEAT  TRANSFER  SUBROUTINES 


The  main  program  of  the  heat  transfer  code  is  a rather  complex  driver 
calling  input  subroutines,  source  generating  routines,  and,  under  their  con- 
trols, the  main  heat  transfer  routine  GETEMP.  Finally,  edit  and  output  rou- 
tines are  called.  The  description  of  that  operation  will  be  given  in 
Section  VII. 

The  present  section  is  devoted  to  the  description  of  the  main  Monte  Carlo 
routine  GETEMP,  which  delivers  a single  estimate  of  the  temperature  at  a given 
point  X in  region  IR,  implementing  the  algorithm  described  in  Section  II. 3, 
and  to  the  descriptions  of  the  main  subroutines  called  by  GETEMP,  which  perform 
the  calculations  described  in  Section  III  and  IV. 

1.  Subroutine  GETEMP 

The  subroutine  GETEMP  is  the  main  Monte  Carlo  routine  which  follows  the 
histories  along  the  lines  described  in  Section  II. 3.  It  delivers  a single  esti- 
mate TEMP  of  the  temperature  at  a point  X at  time  TB  in  region  IR. 

We  first  describe  its  operation  under  normal  circumstances,  which  is  invoked 
when  all  DFAST(IR)  are  set  equal  to  zero.  The  weight  W is  set  to  1.  Subroutine 
GBOX  (described  in  Section  IV. 4.1)  is  called.  A box  is  therefore  constructed 
which  is  wholly  in  region  IR,  and  which  is  centered  at  X.  A call  to  subroutine 
DIFFUS  (described  below)  samples  the  RPP  Green's  function  with  zero  boundary 
conditions  along  the  lines  of  Section  III. 2.  A subsequent  call  to  subroutine 
MOVE  (described  below)  moves  the  point  X to  the  sampled  location  and  updates 
the  time  variable.  If  the  time  variable  is  zero,  the  temperature  at  t=0  in  IR 
is  scored  and  a return  is  executed.  If  t>0,  a test  is  made  whether  the  new  point 
is  on  an  abutting  boundary.  If  it  is  not,  a new  box  is  constructed  around  the 
new  point.  If  this  construction  indicates  that  the  point  was  within  a distance  t 
(defined  in  Section  IV)  of  the  abutting  surface,  the  point  is  treated  as  being 
on  the  abutting  surface.  If  not  on  the  abutting  surface,  the  sequence  of  calls 
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to  DIFFUS  and  MOVE  continues  until  either  the  time  variable  reaches  zero,  or 


I 


the  point  has  diffused  to  a region  boundary. 

If  the  point  is  on  a boundary,  the  neighbor  region  number  IRP  is  deter- 
mined. The  diffusion  coefficient  DP  of  the  neighbor  region  is  looked  up.  The 
convention  is  that  if  DP=0,  the  neighboring  region  is  one  with  known  temperature. 
If  DP<0  then  inhomogeneous  boundary  conditions  apply  with  8=-DP/K  (6  is  defined 
in  Section  II. 2).  If  DP>0,  the  neighboring  region  is  a conductive  one. 

If  DP=0,  a score  is  made,  and  a return  is  executed.  If  DP^O,  a "half  box" 
is  constructed  in  region  IR:  this  is  a box  one  side  of  which  is  on  the  surface 

of  interest,  centered  around  X.  The  box  is  fully  contained  in  IR.  A "half  box" 
is  also  constructed  in  IRP. 

If  DP>0,  the  subroutine  MENDEL  is  called.  Subroutine  MENDEL  is  described 
below;  it  samples  the  two  region  RPP  Green's  function  along  the  lines  discussed 
in  Section  V.  A call  to  subroutine  MOVE  moves  the  point  X and  updates  the  time  t. 
If  t=0,  a score  is  made  and  a return  is  executed.  If  not,  a new  diffusion  loop 
is  started  as  described  above. 

If  DP<0,  the  subroutine  NORMAN  (described  below)  is  called.  It  samples  the 

RPP  Green's  function  with  inhomogenous  boundary  conditions  along  the  lines  of 

Section  III. 5 and  updates  the  weight  value. 

Prior  to  the  call  to  NORMAN,  y is  set,  as  well  as  are  the  estimates  T and 

9 

% 

T needed  for  importance  sampling.  These  are  taken  to  be  the  value  of  the  gas 


temperature  in  IRP  and  the  value  of  the  temperature  in  IR  at  t=0.  If  the 
particle  diffuses  to  the  abutting  boundary,  a score  is  made.  The  same  is  done 
if  t becomes  zero.  In  both  cases,  a return  is  executed.  If  neither  occur,  a 
new  diffusion  loop  is  initiated. 
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The  number  of  diffusion  steps  within  a region  is  limited  to  be  less  or 
equal  to  NRSTEP.  The  number  of  region  crossings  is  limited  by  NRCR.  If  any 
of  these  is  reached,  a return  with  TEMP=-1.  is  executed.  The  same  error  return 
occurs  if  any  of  the  three  dimensions  of  any  box  becomes  smaller  than  DTINY, 
presently  set  to  10  10 . Experience  in  running  the  code  should  teach  how  to 
avoid  such  collapsing  boxes  - a rarely  occurring  event  anyway.  NRSTEP,  NRCR 
are  currently  set  to  200  and  5000,  respectively. 

The  above  description  is  applicable  if  the  input  parameters  DFAST(IR)  is 
set  to  zero  for  all  regions  IR.  It  turns  out  that,  in  solving  heat  transfer 
problems  with  expected  high  temperature  gradients,  homogeneous  regions  have  to 
be  subdivided  into  rather  small  subregions,  to  satisfy  the  requirement  that 
the  diffusion  properties  be  constant  in  each  subregion.  With  such  a detailed 
geometrical  subdivision,  the  calculations  involved  become  excessively  slow, 
because  each  boundary  crossing  involves  two  calls  to  subroutine  GBOX  and  a call 
to  the  boundary  crossing  subroutine  MENDEL.  The  accuracy  of  the  calculations 
involved  in  such  a detailed  treatment  for  crossing  imaginary  boundaries  cannot 
be  justified,  and  an  approximation  to  speed  up  the  calculation  is  in  order. 
Selected  regions  can  be  declared  on  input  as  not  requiring  completely  accurate 
boundary  treatment.  This  is  specified  by  setting  the  input  quantity  DFAST(IR) 
for  selected  regions  IR.  The  approximation  consists  in  assuming  that  the  diffu- 
sion parameters  are  constant  within  a (small)  distance  j DFAST (IR) | from  any 
point  in  region  IR,  including  points  on  the  boundary  of  IR. 

Before  the  call  to  GBOX,  subroutine  GETEMP  examines  the  value  of  DFAST (IR). 

If  zero,  the  calculations  proceed  as  described  above.  If  the  value  of  DFAST (IR) 
is  positive,  the  call  to  GBOX  is  performed  and  the  smallest  dimension  of  the 
RPP  is  compared  to  the  value  of  DFAST (IR).  If  smaller,  the  calculations  are 
performed  as  described  above.  If  larger,  a cube  of  dimension  DFAST (IR)  is  set 
up  and  the  calls  to  DIFFUS,  MOVE,  are  replaced  by  a call  to  subroutine  CUBE. 

If  the  value  of  DFAST (IR)  is  negative,  the  call  to  GBOX  is  bypassed  altogether. 
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a cube  of  dimension  |DFAST(IR) | is  set  up,  and  the  calls  to  DIFFUS,  MOVE, 
are  replaced  by  a call  to  subroutine  CUBE.  Subroutine  CUBE  (described  below) 
moves  the  point  X to  its  new  position,  updates  t,  and,  if  necessary,  changes 
the  region  number  IR. 

2.  Subroutine  DIFFUS 

Subroutine  DIFFUS  samples  homogeneous  RPP  Green's  functions  with  zero 
boundary  conditions,  as  discussed  in  Section  II I. 2. 

2 

The  input  consists  of  a time  cutoff  value  T,  and  of  FOX(I)=a^  /D,  where 
a^,  i=l,3  are  the  dimensions  of  the  RPP,  and  D is  the  diffusion  coefficient. 
It  is  assumed  that  a^  = = 2*D1  and  a^  = 2*D2. 

The  output  consists  of  a time  variable  TMIN^T  and  of  a local  vector 
XX ( I ) <a^/2 . 

The  time  selection  is  achieved  by  three  calls  to  subroutine  PICKT.  The 
position  selection  is  achieved  by  calls  to  subroutine  PICKX. 

3.  Subroutine  MOVE 

This  subroutine  displaces  the  absolute  position  X by  the  displacement  XX 
generated  by  subroutine  DIFFUS.  The  local  vector  XX  is  given  in  a coordinate 

-V 

system  in  which  XX (3)  is  along  the  absolute  vector  Q;  the  remaining  two  co- 
ordinates are  orthonormal  but  otherwise  arbitrary. 

The  time  T is  set  to  T-TMIN. 

4.  Subroutine  CUBE 

Subroutine  CUBE  is  a speeded  up  version  of  the  combination  of  DIFFUS  and 
MOVE,  under  the  simplifying  conditions  that 
FOXX=FOX ( 1 ) =FOX ( 2 ) =FOX ( 3 ) 

and  that  the  vector  Q points  along  the  absolute  z-axis,  i.e.,  that  the  local 
coordinate  system  XX  is  parallel  to  the  absolute  system  X. 
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5.  Subroutine  NORMAN 

Subroutine  NORMAN  is  the  equivalent  of  DIFFUS,  except  that  homogeneous 
boundary  conditions  are  imposed  at  XX(3)=0.  Importance  sampling  is  implemented 
and  the  calculations  follow  the  lines  described  in  Section  III. 7.  The  selection 
of  times  t^  and  t2  is  achieved  by  calls  to  subroutine  PICKT.  The  selection  of 
t3  is  via  a call  to  a subroutine  PICKTG.  Similarly,  selection  of  XX (1)  and 
XX (2)  is  via  calls  to  PICKX,  whereas  the  selection  of  XX (3)  i3  via  a call  to  a 
subroutine  PICKXG. 

6.  Subroutine  MENDEL 

This  is  the  boundary  crossing  subroutine.  First  the  diffusion  coefficients 
and  the  heat  conductivities  in  the  two  regions  are  compared.  If  both  comparisons 
show  agreement  within  5%,  an  average  diffusion  coefficient  is  calculated,  and 
a call  to  DIFFUS  is  followed  by  a return. 

If  the  homogeneity  test  fails,  the  subroutine  switches  to  the  method  de- 
scribed in  Section  IV.  The  part  of  the  calculation  involving  the  floating  RPP 
method  is  performed  in  subroutine  DIMOR. 

7.  Subroutine  DIMOR 

Constructs  an  RPP  within  an  RPP,  and  samples  time  and  position  via  calls 
to  PICKT  and  PICKX. 

8.  Auxiliary  Subroutines 

Special  subroutines  perform  sampling  of  one  dimensional  Green's  functions, 
as  described  in  Section  III.  These  include: 

PICKT,  PICKX  - See  Section  III. 4 
PICKTG,  PICKXG  - See  Section  III. 6. 2. 

When  appropriate,  the  early  time  approximation  is  invoked  by  calls  to  sub- 
routines EARLTG,  EARIiXG.  (See  Section  III. 6.1.)  When  appropriate,  the  eigen- 
values are  calculated  by  a call  to  subroutine  TRANS.  The  function  X=GAUSS(I) 
returns  X>0  sampled  from  2/JH  exp(-x2)dx.  The  function  y=ECHERF(X)  returns 
y=2//rr  exp(X2)  JX  exp(-y2)dy. 
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VII.  THE  THREE  DIMENSIONAL  TIME  DEPENDENT 


ADJOINT  MONTE  CARLO  HEAT  TRANSFER  CODE 

The  simplest  application  of  the  code  is  to  a problem  where  the  diffusion 
properties  of  the  medium  are  independent  of  the  local  temperature.  In  that 
case,  the  operations  of  the  code  are  reduced  to  reading  the  geometrical  descrip- 
tions, the  initial  conditions,  the  boundary  conditions,  the  coordinates  of  de- 
tector points,  the  value  of  the  time  variable  at  which  the  temperature  needs 
to  be  calculated,  and  the  number  of  Monte  Carlo  histories  to  be  followed.  The 
code  then  performs  that  many  calls  to  subroutine  GETEMP,  calculates  the  effected 
temperature  and  its  variance,  and  prints  the  results. 

A more  complex  problem  arises  in  the  case  where  the  diffusion  properties 
do  depend  on  the  local  temperature.  The  time  span  between  initial  and  final 
time  must  then  be  subdivided  into  a number  of  time  bins.  The  time  bins  must  be 
small  enough  so  that  the  local  diffusion  properties  can  be  considered  as  constant 
throughout  the  time  bins.  Homogeneous  regions  must  be  subdivided  into  subregions. 
Each  subregion  must  be  small  enough  so  that  the  local  temperature  at  any  point 
in  the  subregion  does  not  differ  appreciably  from  the  average  temperature  in 
the  region. 

Given  such  a subdivision  of  space  and  time,  the  operation  of  the  code  is 
as  follows.  The  initial  conditions  are  read  for  the  very  first  time  bin  only. 

The  boundary  conditions  are  also  read  in.  The  diffusion  properties  are  calcu- 
lated for  each  region  at  the  value  of  the  initial  temperatures.  These  are  assumed 
to  remain  constant  throughout  the  time  bin.  The  code  then  estimates  temperatures 
at  the  end  of  the  time  bin  for  all  regions.  These  calculated  temperatures  serve 
as  initial  conditions  for  the  next  time  bin. 

The  calculation  of  region  temperatures  is  performed  in  the  following  order. 
The  code  searches  for  the  region  with  the  highest  known  temperature  at  the  end 
of  the  time  bin.  (This  is  restricted  to  external  regions  only  at  the  beginning 
of  the  calculation.)  Its  neighbors  are  examined.  The  first  neighbor  found  for 
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which  the  temper ature  has  not  been  calculated  yet  is  then  picked  up,  and  the 

temperature  is  estimated  by  adjoint  Monte  Carlo.  When  all  the  neighbors  have 

been  completed,  the  region  is  excluded  from  the  list,  and  the  next  region  with 

the  highest  known  temperature  is  searched.  The  calculations  continue  as  long 

as  such  highest  temperature  is  appreciably  higher  than  T , where  T is  set 

m m 

slightly  higher  than  the  lowest  temperature  of  the  configuration  at  the  be- 
ginning of  the  time  bin.  All  remaining  regions  are  surrounded  by  regions  with 

temperature  less  than  T ; their  temperature  is  set  equal  to  T . 

m z 

The  calculation  of  temperatures  averaged  over  regions  is  performed  by  cal- 
culating the  temperature  at  points  distributed  uniformly  over  the  region;  the 
sampling  of  such  points  is  described  in  Section  VII. 5 below.  The  efficiency  of 
that  sampling  may  be  improved  by  providing  additional  input  as  described  there. 

1.  The  Main  Program  HEATON 

The  main  program  is  a driver  which  reads  input,  and,  under  its  control, 
calls  additional  input  and  processing  routines,  and  calculational  routines. 

A complete  description  of  the  input  i3  given  in  Section  IX.  The  first 
card  read  provides  the  main  specifications:  NHIST  is  the  number  of  histories  to 
be  run  for  each  region,  for  each  time  bin  (or  for  each  point  if  point  detectors 
are  specified) . NSTAT  (<200)  is  the  size  of  a group  of  histories  for  the  pur- 
pose of  source  generation.  ITCUT  is  a CPU  time  cutoff:  the  calculations  will  be 
properly  terminated  at  the  end  of  a time  bin,  and  a restart  tape  will  be  written 
if  CPU  time  is  nearing  that  limit.  IT1  is  the  first  time  bin  to  be  considered 
(IT1=1  at  first,  if  no  restart  tape  is  available).  IT2  is  the  last  time  bin 
to  be  considered.  If  IPD*0,  the  last  time  bin  IT2  is  treated  like  any  other  ones. 
If  IPD>0,  the  temperatures  are  calculated  only  at  the  given  IPD  point  detectors 
for  time  bin  IT2.  IREZ  should  normally  be  equal  to  zero.  If  set  equal  to  1, 
geometry  input  will  not  be  expected,  but  geometry  information  will  be  picked  up 
from  a restart  tape  even  if  ITl=l.  IPLO  is  a flag  controlling  the  quantity  of 
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printed  output  and  calls  to  subroutine  PLOT.  If  IPLO  is  even  (including 
zero  or  blank),  all  input  quantities  are  printed  back.  If  odd,  the  bulk  of  the 
output  is  suppressed,  except  for  estimated  temperatures.  Subroutine  PLOT  will 
be  called  if  IPLO  is  equal  to  2 or  3.  LRN  controls  the  initialization  of  the 
random  number  routine.  Leaving  that  as  blank  will  leave  the  random  number  as  1 
if  IT1=1 , or  as  recorded  on  the  restart  tape.  If  non-zero,  it  must  be  an  odd 
integer. 

The  geometry  input  (if  called  for)  is  read  in  by  subroutine  GENI  (see 
Section  V) . The  coordinates  X,  Y,  Z,  of  a point  within  the  configuration  must 
also  be  provided. 

The  next  item  of  input  deals  with  reflection  regions  and  translationally 
repetitive  arrays.  These  are  described  in  Section  VII. 9. 

Optional  input  to  improve  the  efficiency  of  source  generation  is  read  in  by 
the  main  program.  The  meaning  of  that  input  will  become  apparent  in  Section  5. 
For  those  regions  for  which  no  input  has  been  specified,  or  for  which  none  is 
available  from  the  restart  tape,  the  information  is  generated  by  subroutine 
SOUSET.  If  either  input  or  generated,  the  information  is  checked  by  subroutine 
SOUCK. 

The  physical  compositions  are  defined  next.  The  temperature  dependence  of 
the  diffusion  parameters  is  defined  as  follows: 

Product  of  density  and  specific  heat: 

2 

PC  = rQ  + r2  T 

Thermal  conductivity 

K = kQ  + k2  T2,  where  k^O. 

The  input  consists  of  the  parameters  rQ,  r2,  kQ,  k2  for  the  NCOMP  materials  pre- 
sent in  the  problem.  The  first  material  entered  is  assumed  to  be  the  main  one. 
The  transformations  discussed  in  Section  II. 5 are  performed  on  its  basis. 

The  meaning  of  the  quantities  DFAST  is  described  in  Section  VI. 1. 


r,  • 

L’ 


> 


If  IT1«1,  the  initial  conditions  are  read  in.  The  reading  of  the  boundary 
conditions  is  done  by  subroutine  REATIM. 

If  IT2=*0,  the  coordinates  x,y,z  of  point  detector  are  read  in  and  the  sub- 
routine REGCAL  is  called  to  calculate  the  temperature  at  that  point. 

If  IT2>ITl,  the  calculation  of  region-averaged  temperatures  is  done  by 
subroutine  TIMSTP . 

Two  restart  tapes  are  generated;  TAPE  8 involves  processed  geometry  informa- 
tion, including  information  learned  about  neighbors,  and  source  generation  in- 
formation. TAPE  9 involves  the  temperatures  at  the  end  of  the  time  step. 

2.  Subroutine  REATIM 

The  subroutine  reads  the  parameters  of  boundary  conditions  at  the  end  of 
the  current  time  step.  Linear  interpolation  is  assumed  between  the  beginning  and 
the  end  of  the  time  step. 

3.  Subroutine  TIMSTP 

The  subroutine  searches  for  the  region  with  the  highest  known  temperature, 
examines  its  neighbors,  and,  for  each  region  for  which  a calculation  is  needed, 
calls  the  subroutine  REGCAL. 

4.  Subroutine  REGCAL 

The  subroutine  either  calls  the  source  generating  routine  SOUSET,  or  picks 
up  the  coordinates  of  a detector  point,  and  performs  NHIST  calls  to  subroutine 
GETEMP  (see  Section  VI. 1).  It  calculates  the  mean  and  root  mean  square  tempera- 
ture and  prints  one  line: 

JR  = region  number 

TQ  = mean  temperature  T 

TR  = root  mean  square  temperature 

TR  = mean  0 

TS  = root  mean  square  0 

CPU  = CPU  time  in  seconds  for  that  calculation. 


see  Section  II. 5 
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The  subroutine  produces  NSTAT  points  uniformly  distributed  in  region  IR, 
under  the  asusmption  that  the  region  is  completely  enclosed  by  the  sector  of 
a circular  cylinder  defined  by  the  two  figures  above.  The  relevant  parameters 


XSTART (1 / IR)  = x. 


XSTART ( 2 , IR) 


XSTART (3 ,IR)  = z 


central  point  on  axis 


XSTART (4, IR)  = R Outer  radius 
XSTART (5, IR)  = r Inner  radius 


XSTART (6, IR)  = W Half  Height 


AZIM(1,IR)  **  Hi  = Minimum  azimuth 

AZIM (2 , IR)  = 1^2  = Maximum  azimuth 
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The  code  samples  points  uniformly  distributed  in  the  sector,  and  rejects 
the  sample  if  the  point  is  outside  of  IR.  In  the  process,  the  subroutine 
checks  for  point  in  IR  with  a distance  to  the  axis  which  is  either  greater  than 
R or  smaller  than  r.  If  such  points  are  found,  R and  r are  properly  updated, 
and  points  previously  generated  within  the  current  aggregate  are  properly 
corrected  if  necessary.  No  such  check  is  performed  on  W,  i p . 

The  efficiency  of  the  routine  depends  on  the  fit  of  the  sector  to  region 

IR. 

6.  Subroutine  REGSET 

Subroutine  REGSET  generates  the  XSTART  parameters  for  those  regions  for 
which  that  input  has  not  been  provided.  The  minimum  azimuth  is  set  equal  to  0, 
and  the  maximum  azimuth  to  2ir;  thus  only  full  cylindrical  annuli  are  generated. 

The  first  region  considered  is  one  for  which  any  internal  point  P=(x,y,z) 
is  known.  If  none  other,  it  may  be  the  point  read  in  by  the  main  program  after 
the  subroutine  GENI  has  been  called.  A ray  is  fired  isotropically  from  P and 
points  of  intersection  of  this  ray  with  region  IR  are  generated.  The  point  P 
is  moved  to  the  midpoint  between  the  pair  of  points  which  are  most  distant  from 
each  other.  A new  ray  is  fired  from  the  new  P position,  and  intersections  are 
added  to  the  list....  The  process  terminates  when  200  points  have  been  gener- 
ated on  the  surface  of  IR.  This  population  of  points  is  then  fitted  by  the 
"best"  cylindrical  annulus  in  the  z-direction. 

The  fitting  is  actually  to  the  smallest  finite  cylinder  in  the  z-direction. 
The  cylinder  extends  from  the  smallest  z-coordinate  to  the  largest  z-coordinate. 
The  smallest  circle  enclosing  the  population  of  x,  y points  is  determined.  This 
is  done  as  follows.  A first  pass  through  all  pairs  of  points  determines  the 
pair  of  points  x^,x2  with  the  largest  separation  distance.  A tentative  circle 
of  radius  equal  to  one  half  of  that  distance  is  centered  at  the  midpoint  of 
that  pair.  If  all  points  are  within  that  circle,  the  tentative  circle  is  the 


smallest  possible  circle-  To  find  out,  a next  pass  through  all  the  points  de- 
termines the  point  x^  with  the  largest  distance  to  the  tentative  center.  If 
that  distance  is  larger  than  the  tentative  radius,  the  circle  is  drawn  as 
passing  through  the  points  x^,  x2,  x^.  Although  this  last  circle  is  not  neces- 
sarily the  smallest  possible  one,  it  is  kept  as  is. 

During  the  process  of  ray  firing,  the  neighbor  IRP  of  region  IR  is  deter- 
mined at  each  intersection.  If  no  starting  point  has  yet  been  determined  for 
IRP,  the  intersection  point  (slightly  displaced)  is  stored  as  one. 

7.  Subroutine  SOUCK 

This  subroutine  performs  the  same  operations  as  subroutine  REGSET,  but 
with  many  more  checks.  It  also  keeps  track  of  the  efficiency.  The  content  of 
the  XSTART  and  AZIM  arrays  is  printed.  If  the  efficiency  becomes  intolerable, 
the  holler ith  information  BAD  is  also  printed,  thus  suggesting  that  source  in- 
formation be  specified  on  input. 

8.  Subroutine  PLOT 

This  subroutine  provides  a rather  limited  capability  of  obtaining  printer 
plots  of  x-y  cuts  thro  igh  the  geometry. 

9.  Reflection  Regions  and  Translationally  Repetitive  Arrays 

Only  regions  described  as  single  planes  can  be  declared  as  reflection 
regions.  If  the  corresponding  composition  number  ICOMP  is  set  to  any  positive 
number,  reflection  boundary  conditions  will  be  imposed  on  the  plane  defining 
that  region. 

The  geometrical  description  of  a configuration  consisting  of  groups  of 
regions  called  cells  which  are  translationally  repetitive  can  be  simplified  to 
the  description  of  only  one  cell  with  the  use  of  translationally  repetitive 
arrays.  The  conditions  for  repetitive  geometry  are  given  by  the  following. 
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Given  a vector  D,  a repetition  number  m and  a total  number  of  regions  n*m, 
if  a point  X is  in  region  i then  point  (X+S)  is  in  region  (i+n) , provided 
i+n<n-m,  and  point  (X-D)  is  in  region  (i-n) , provided  i-n>0  (see  Figure  1). 

For  the  configuration  shown  in  Figure  2,  m repetitive  cells  can  be  de- 
fined. To  describe  the  entire  configuration,  only  the  primary  cell  needs  to 
be  defined.  Thus,  in  Figure  3,  consider  the  first  cell  which  can  be  described 
by  n regions.  The  boundary  surfaces  are  described  by  regions  i^  and  i^,  such 
that  0<i^<n,0<i2<n.  The  rest  of  the  cell  can  be  arbitrarily  subdivided  into 
(n-2)  regions  numbered  1,  2,  ...n  excluding  i^  and  i2.  The  configuration  is 
composed  of  m of  these  cells  so  that  the  external  configuration  surfaces  are 
defined  by  regions  i^  and  i2  + (m-1) -n.  The  boundary  conditions  at  these  sur- 
faces are  determined  by  the  signs  of  JCOMPfi^)  and  ICOMP (i2+ (n-1) *n) . If 
ICOMP (J)<0,  the  usual  conventions  hold,  i.e.,  O-known  temperature,  <0-radia- 
tion  types;  if  ICOMP (I)>0,  reflection  boundary  conditions  are  imposed. 

The  input  to  describe  repeating  arrays  and  reflection  regions  consists 
of  a list  of  up  to  10  region  numbers.  The  first  two  regions  correspond  to 
regions  i^  and  i2  defined  above,  i.e.,  regions  on  either  side  of  the  first 
cell.  The  remaining  regions  in  the  list  are  reflection  regions  in  the  first 
cell.  The  repetition  number  is  also  read  in.  If  there  is  only  one  region 
number  in  the  list  or  if  m=l,  then  every  region  in  the  list  is  a reflection 
region. 

With  the  use  of  repeating  arrays,  the  region  and  helps  (ITEM  6)  input  is 
reduced  in  that  only  one  cell  of  n regions  need  be  described.  Computer  memory 
requirements  are  also  reduced.  Temperature,  composition  and  DFAST  array  input 
still  consists  of  data  for  all  n.m  regions. 
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TABLE  2 


DIMENSION  OF  LARGE  ARRAYS 


COMMON  BLOCK 

ARRAY  NAME 

DIMENSION 

BLANK 

I BLOC 

NDQ1 

ISLOC 

NDQ1 

MAS 

NDQ 

SOU 

xs 

3*ni 

XSTART 

6»NDQ2 

AZIM 

2-NDQ2 

INP 

DIFF 

NDQ  3 

AAK 

NDQ3 

TO 

NDQ3 

T1 

NDQ3 

RES 

I COMP 

NDQ  3 

FAST 

DFAST 

NDQ  3 

DI 

ALP 

n2 

ALS 

n2 

GG 

n2 

CC 

n2 

GQ 

3 "n2 

NDQ  = large  enough  to  store  geometry  information  and  neighbor  information 
NDQ1  NS  = number  of  surfaces 
NB  = number  of  bodies 
NR  = number  of  geometrical  regions 
NDQ 2 ^ NR  = number  of  geometrical  regions 
NDQ3  NPHYS  = NBI*NR  = number  of  physical  regions 
n^  = 200  = maximum  size  of  statistical  aggregates 
n^  = 200  = maximum  number  of  eigenfunctions 


— 


IX.  INPUT  DESCRIPTION 


The  card  input  consists  of  the  following  items: 

Item  1 FORMAT  (715,120) 

NHIST  Number  of  histories 

NSTAT  Group  of  histories  (<200) 

ITCUT  Real  time  cutoff 

IT1  First  time  bin  to  be  treated  (>1) 

IT2  Last  time  bin  to  be  treated  (>ITl) 

IPD  Number  of  point  detectors 

IREZ  Restart  option  switch  (1  on,  0 off) 

IPLO  Output  control:  input  printed  back  if  IPLO=0  or  2 

plot  routine  called  if  IPLO=2  or  3 

LRN  Random  number  initialization  (0-use  default  - must  be 

odd  positive  otherwise) 

Omit  Items  2-4  if  IT1>1  or  IREZ^O 


Item  2 GEOMETRY  INFORMATION 

Item  2.1  Header  card  FORMAT  (I4,19A4) 

IPR  If  zero,  input  and  processed  data  printed 
HOLL  Any  hollerith  information 
Item  2.2  Surface  input  FORMAT  (2X,A3,8E9. 3) 

TYP  Surface  type  identifier 


EPS 

X 

Y 

Z 


"Surface  thickness" 
coordinates  of  vertex 


A^-A^  - as  specified  in  the  following  table: 


‘..-r 

k ■< I 
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Item  2.3  Body  input 


FORMAT  (2X, A3, 1515) 


TYP  Body  type:  only  GEN  allowed.  END  is  end  of  list. 


List  of  signed  surface  numbers 


Item  2.4  Region  input  FORMAT  (2X, A3, 1515) 


ANY  Any  three  characters.  END  is  end  of  list. 


List  of  signed  body  numbers 


Item  3 FORMAT  (3E15.5) 


Y > Coordinates  of  any  point  within  the  configuration 

2 J 

Item  4 Repeating  Array  - Reflector  regions 
Item  4.1  FORMAT  (1015) 

NREF  No.  of  regions  given  (NREF<10) 


Region  numbers 


NBI  Repetition  number 


Note  - Omit  Item  5 if  IPLO=0  or  1 
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Item  5 x-y  plot  input 


FORMAT  (5E12.4) 


For  each  plot  enter 


coordinate  of  left-upper  corner 


2 f coordinates  of  right  lower  corner  (x^x^,  y 2<y^) 


A blank  card  terminates  plot  input 


Item  6 Optional  input  for  source  generation 
Item  6.1  FORMAT  (15) 


Number  of  regions  for  which  input  is  provided 


Note  - omit  Item  6.2  if  NHE=>0 


Item  6.2  FORMAT  (I5,10E7.0) 


For  NHE  regions,  enter  a card  consisting  of  the  description: 

IR  Region  number 

The  following  defines  a sector  of  a cylindrical  annulus  parallel  to  the 
z-axis,  which  encloses  region  IR: 


Coordinates  of  midpoint  on  axis 


R Outer  radius  of  cylinder 
r Inner  radius  of  cylinder 
w Half  height  of  cylinder 


A cos^j 
A c ost)> 


A arbitrary  - ■ minimum  angle 

If  A=0,  is  set  to  0 


B cosily  I B arbitrary  - t|»2  » maximum  angle 

B sin^2  J If  B=0,  ip 2 is  set  to  2ir 

(Notation  refers  to  Figure  in  Section  VII. 5) 


* Information  for  selected  regions  can  be  entered  in  any  order.  The  information 
currently  entered  for  region  IR  will  be  recorded  on  a restart  tape.  That  in- 
formation will  overide  information  recorded  on  a previous  restart  tape  (if  any), 


Item  7 


COMPOSITION  INFORMATION 


Item  7.1  FORMAT  (15) 


NCOMP:  Number  of  compositions 


Item  7.2  FORMAT  (4E15.5) 


For  each  of  the  NCOMP  compositions,  enter: 


parameters  defining  the  product  of  density 

2 

and  heat  capacity  pC  = RCQ  + RC2  ? T 
parameters  defining  the  heat  conductivity 


K = K0  + K2 


(The  first  composition  entered  is  the  main  composition  of  the  problem. ) 

Item  7.3  FORMAT  (1015) 

(ICOMP (IR) ,IR=1,NR)  For  regions  internal  to  the  configuration,  enter 

composition  number.  For  external  regions,  enter  0 if 


temperature  is  known,  and  -1  if  boundary  condition  of 
the  radiation  type  applies.  Enter  any  positive  number 


for  reflection  regions. 

Item  8 INFORMATION  FOR  "SMALL  BOX  APPROXIMATION"  FORMAT  (3E10.2) 

(DFAST (IR) ,IR=1,NR)  Enter  zeroes  for  no  approximation.  See  text,  otherwise: 


if  DFAST <0 , diffuse  with  box  of  fixed  size  |DFAST|  . If 
DFAST>0  diffuse  with  either  actual  box,  or  box  of  size 
DFAST,  whichever  is  larger. 


Note  - Omit  Item  9 if  IT1>1 


Item  9 INITIAL  CONDITIONS 


Item  9.1  FORMAT  (E15.5) 


TBl : Initial  value  of  the  time  variable 


Item  9.2  FORMAT  (5E15.5) 


(Tl (IR) ,IR=1,NR) : Initial  temperatures  for  all  regions. 
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Note  - Repeat  items  10-12  for  each  time  bin 


1 


: 


l>  1 

1 ' 
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Iteu.  10  FORMAT  (E15.5) 


TB:  value  of  the  time  variable  at  end  of  time  bin 

Item  11  Boundary  conditions  FORMAT  (5E15.5) 

T2 ( IR1 ) Final  temperature  of  first  region  with  ICOMP (IRl)^O 

T2 ( IR2 ) Final  temperature  of  second  region  with  ICOMP (IR2)^0 

T2 ( IRL)  Final  temperature  of  last  region  with  ICOMP (IRL)^O 

Item  12  Coefficients  of  surface  heat  transfer.  FORMAT  (5E15.5) 


H(IR1) 
H (IR2) 


Enter  value  of  coefficient  if  ICOMP  (IRN) <0. 
Leave  blank  if  ICOMP (IRN)=0 


H (IRL) 

Note  - Omit  Item  13  if  IPD=0 

Enter  Item  13  IPD  times 
Item  13  Detector  input  FORMAT  (3E15.5) 


x 

y 

z 


] 


Coordinates  of  the  point  detector 


Item  14  END  OF  FILE 


Input  Tapes 

Tape  5 - Card  input 

Tape  8 - Geometry  restart  tape  - needed  if  IT1>1,  or  if  ITl=l  and  IREZ=1. 
Tape  9 - Initial  conditions  - needed  if  ITl>l. 

Output  Tapes 

Tape  6 - Printer 

Tape  8 - Geometry  restart  tape 

Tape  9 - Final  conditions  - generated  if  IT2>1  or  if  IT2=1  and  IPD=0. 
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